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Chapter  1 

Introduction 


The  BTEC  thermal  model  is  a  1-D  and  2-D  cylindrical  coordinate  system  simulation  of  optical  radiation  and  radio 
frequency  thermal  interaction  with  tissues.  The  code  supports  the  illumination  of  tissues  by  sources,  the  temperature 
response  from  the  linear  absorption  of  optical  radiation,  and  analysis  of  subsequent  damage.  The  model  takes  its  name 
from  the  four  primary  authors  of  the  model:  Buffington,  Thomas,  Edwards,  and  Clark,  although  many  others  have 
contributed. 

Finite  difference  numerical  methods  are  employed  in  the  solutions  of  heat  transfer.  An  alternating-direction-implicit 
(ADI)  finite  difference  method  is  employed  in  the  solution  of  the  heat  equation  in  2-D,  while  a  Crank-Nicholson 
Method  is  employed  in  the  1-D  solution.  The  BTEC  simulation  can  be  configured  with  a  source  term  defined  by  a 
single  or  by  multiple  emitters,  such  as  laser,  RE,  or  broadband.  The  BTEC  can  also  compute  a  source  based  on  the 
initial  electric  field  profile  and  a  refractive  index  distribution. 

The  tissue  simulation  is  represented  as  a  one-dimensional  stack  of  homogeneous  layers  along  the  z-axis.  Each  layer 
can  have  differing  thermal,  optical,  and  physical  properties.  The  linear  absorption  coefficient  of  each  layer  defines  the 
energy  transfer  from  the  optical  source  to  the  tissue.  Boundary  conditions  along  the  axial  and  radial  (2-D)  coordinates 
may  be  selected  as  a  sink  (temperature  held  constant)  or  as  a  combination  of  surface  boundary  conditions  (convection, 
radiative,  or  evaporative). 

The  BTEC  model  currently  employs  a  single  rate-process  model  of  thermal  injury.  User-defined  parameters  for  dam¬ 
age  rates  and  activation  energies  as  a  function  of  temperature  can  be  programmed  for  each  tissue  type  or  layer  used.  A 
number  of  damage  integral  values  or  temperature  shift  searches  are  available  in  order  to  estimate  damage  thresholds 
or  for  comparison  to  experimental  data. 
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Chapter  2 


One-Dimensional  Heat  Equation 


2.1  Heat  Equation 

The  BTEC  model  employs  a  Crank-Nicholson  Method  [9]  for  the  solution  of  the  time-dependent,  one-dimensional 
heat  transfer  problem.  A  detailed  derivation  of  the  method  is  presented  below. 


The  formulation  begins  with  the  time-dependent  heat  equation  expressed  in  equation  2.1: 

dv(r  t ) 

p(r)  c(r) — =  V  ■  Mr,  v)Vv(r,  f)]  +  g(r,  t )  (2.1) 

Here,  p( r)  is  the  material  density,  c(r)  is  the  specific  heat,  and  k( r,  v)  is  the  thermal  conductivity.  The  expression  v(r,  t) 
is  the  temperature  rise  and  g( r,  t )  is  the  source,  or  driving  term.  The  source  term  contains  the  effect  of  an  optical  or 
radio-frequency  source  along  with  the  effects  of  perfusion  (blood  flow). 


The  product  rule  is  used  to  expand  this  equation,  yielding 

dv(r,  t )  , 

p( r)  c(r) — - —  =  V/r(r,  v)  •  Vv(r,  t)  +  k( r,  v)  V  v(r,  v)  +  g(r,  t). 
at 

equation  2.2  can  then  be  expanded  for  one-dimensional  Cartesian  coordinates  using  the  following: 

_  Ok  „ 

V/r  =  — x 
dx 

dv  ^ 

Vv  =  — x 
dx 

2  S'-v 

Tv  =  a? 

Utilizing  Equations  2.3,  2.4,  and  2.5,  the  expanded  heat  equation  is 


dv(x,t)  dK(xJ)dv(x,t)  d2v(x,  t) 

p(x)c{x) - - -  =  - - - I -  K(X,t) —5 - \-g{x,t). 


dt 


dx  dx 


dx2 


(2.2) 

(2.3) 

(2.4) 

(2.5) 

(2.6) 
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2.2  Finite  Difference  Representation 


The  expanded  heat  equation  from  equation  2.6  is  transformed  to  a  finite  difference  representation.  Space  is  discretized 
such  that  calculations  are  computed  from  an  indexed  position,  x,  where 

i  =  0, 1. 

This  spatial  representation  yields  a  finite  difference  representation: 


PiCr 


A  t 


v"  x  dxv"+1+6y;  62y;+1+62xy 

-  =  OxKi  - - -  +  Kj - - - 


+  8i 


(2.7) 


In  equation  2.7,  the  Crank-Nicholson  Method  is  used.  The  forward  finite  difference  at  time  tn  and  backward  finite 
difference  at  time  tn+l  are  averaged  in  the  Crank-Nicholson  Method. 


By  multiplying  equation  2.7  by  2  the  expression 


2piCj  /„+ 1 

At  '  '■ 


v”)  =  6xk i  6XVI+1  +  6xKi  6XV1  +  Ki  S2v?+1  +  *,■  Sty  +  2 ,g{ 


(2.8) 


is  found.  Next,  to  group  known  and  unknown  terms,  all  v'!+1  terms  are  moved  to  the  left  hand  side  (LHS)  and  all  v" 
terms  are  moved  to  the  right  hand  side  (RHS): 


2 PiCi  >i+l 

Ar 


(  8xki  6xVl+x  +  Ki82X+l) 


2-pic  j 

At 


Vj  +  dxKj  6XV‘  +  Ki  82y;  +  2 gi 


(2.9) 


The  derivative  operators  on  v  are  expanded  for  final  simplification  using  stretchable  finite  differences.  Stretchable 
finite-derivatives  are  used  if  a  non-uniformly  spaced  grid  is  desired.  Finite  differences  are  used  to  represent  continuous 
derivatives  in  discrete  computational  space.  Table  2.1  shows  finite  difference  operators  that  will  be  used.  Due  to  the 
size  of  the  equations,  the  calculations  are  presented  in  two  parts,  the  LHS  and  the  RHS. 


Operator 

Un-stretched 

Stretched 

Sxfi 

1 

+ 

A  Xj-  r  AX/+-A*;_  r  AjC/+  r 

2  h 

AxjAxj+  Jl+  ^  Axj-Axj+  Jl  AxjAxi-  Jl~^ 

fM-lf+f-i 

2  f  2  r  .  2  r 

h2 

Axi&x^I'+1  Ax^&x^I'  1  AxiAxi-J'-1 

Table  2.1:  Finite  Difference  Operator  Table  for  1-D  Heat  Equation 
Here  h  is  a  constant  step  size  between  points  and  the  Ax’s  are  debited  as 


Axi  =  x,+i  -  X;_  i 


Ax,+  =  x,+i 


-  x, 


Ax^  —  Xj 


Xi-l. 
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Left-Hand  Side  of  1-D  Finite  Difference  Form 

The  LHS  is  defined  as 


Expanding  the  operators  on  v. 


LHS  i  = 


2pici  1 

At  ‘ 


(  SxKi  Sxv"+'  +  Ki62xVI+l)  . 


LHS  i  =  ^v?+1 


dxKi 


At  1 

Ax, 


Ax, Ax, 
2 


^+1  +  Ax,+  -  Ax,_^1+1  _  Ax,+  ^+1 
!+1  Ax,_Ajc,+  1  AxjAx,-  ,_1 


i+ 


,«+i 


- 


Ax,Ax,+  1  Ax;_Ax;+ 

and  collecting  similar  v  terms,  the  expression 


2  _y"+1  +  ^  _„+l 


AxiAxt. 


LHS ,  =  - 


1 


[2k,  +  ^k,Ax,_]  v". 


AxjAx/ 

IpiCi  6xKj[Axi+  -  Ax,_]  -  2 Kj 


:+l 
+  1 


Af 

1 


Ax, Ax, 

is  found.  Now,  if  the  following  constants  are  defined. 


Ax,_Ax,+ 

[2k,  -  <1xk,Ax,+]  vft1 


v? 


+i 


1 


[2k,  -  6xKiAxi+] 


L  2 pid  6xk,[Ax1+  -  Ax,-]  -  2k, ■ 

b:  = 


At 


cf  =  - 


1 


AxjAx, . 


Ax,_Ax,H 


[2k, ■  +  dtK,Ax,_] 


the  LHS  can  be  represented  in  a  simplified  format: 

LffS^afv-Z+^'+cfy;.:/ 


Right-Hand  Side  of  1-D  Finite  Difference  Form 

Following  the  same  formulation  as  for  the  LHS,  the  RHS  becomes 

RHSi  =  ^  +  5xKi  6xv;  +  Ki  62xv;  +  2gi. 

The  operators  on  v  are  expanded  to  form 


(2.10) 


(2.11) 


(2.12) 


(2.13) 

(2.14) 

(2.15) 


(2.16) 


(2.17) 


RHS,  = 


5xk t 


+  Ki 


At 
A  x,-_ 
Ax, A  x,+ 
2 


■  v»  +  Ax,,  -Ax,  ^ 
1+1  Ax,_Ax,+  ' 
2 


Ax,, 


Ax, A  X,-. 


v?+t  - 


Ax,  A  x,-_  '  1 
2 


Axj-Ax; 


-v  + 


+  2gi 


(2.18) 
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and  then  by  collecting  similar  v  terms,  the  expression 


is  found.  Again,  the  constants 


RHSi  = 


1 


Ax,Ax,+ 


[2 Ki  +  8xKiAxi-\  <+i 


2 ptCj  8xKi[Axi+  -  Ax,_]  -  2 ^ 


At 

1 


Ax;-Ax;j 


v? 


AxjAxj- 


[2 Ki  -  6xKiAxi+ \  V\_x  +  2gj 


[2 Ki  -  8xKjAxi+] 


R  _  2 pid  i  SxKt[ Axi+  -  Ax,_]  -  2/q 


b?  = 


At 


Ax/_Ax/+ 


1 


Ax,  A  x/. 


[2k;  +  dxX;Ax,^] 


are  defined  so  that  the  RHS  can  be  simplified  to 

RHSi  =  afv?_i  +  bfv1;  +  c?v?+1  +  2ft. 
Combining  Equations  2.16  and  2.23  the  simplified  set  of  equations  appears  as 


«  +  =  afvtr  +  +  cf  v"+1  +  2g, 


(2.19) 

(2.20) 

(2.21) 

(2.22) 

(2.23) 

(2.24) 


2.3  Solution  Matrix  Representation 


The  RHS  of  equation  2.24  can  be  found  explicitly  everywhere  and  is  a  constant  to  be  calculated.  For  this  derivation 
when  a  boundary  is  encountered,  the  edge  point  will  be  held  at  zero.  This  is  the  simplest  boundary  condition  and  will 
hold  the  edge  values  to  zero.  In  the  future,  these  edge  points  will  need  to  be  modified  for  other  boundary  conditions. 
The  LHS  of  equation  2.16  must  be  solved  implicitly  using  the  Thomas  algorithm  [9],  Below  in  equation  2.25  is  the 
matrix  representation  of  the  tridiagonal  system  that  must  be  solved: 


0  a[ 


0 


0 


1 

0 


bL 

UM-\ 

aM 


''M- 1 

bL 

UM 


c1 

k',+1 

v"+1 

.,n+ 1 


n+l 
VM- 1 
H  +  l 
VM 


RHS  o 
RHS  i 
RHS  2 

RHSi 


RHS  m  i 
RHS  M 


(2.25) 


In  equation  2.25,  it  is  assumed  that  v_i  =  0  and  vm+i  =  0.  This  is  the  simplest  boundary  conditions  corresponding  to 
a  sink.  Other  boundary  condition  are  possible  by  modifying  the  terms  bo,  c o,  bM,  RHS o,  and  RHS  m- 
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2.4  One-Dimensional  Boundary  Conditions 


2.4.1  Constant  Temperature  Boundary  Condition 

Also  known  as  the  sink  boundary  condition,  this  is  the  simplest  of  the  boundary  conditions.  This  boundary  condition 
holds  a  boundary  temperature  change  equal  to  zero.  In  this  model,  the  constant  temperature  boundary  condition 
implies  that 


v 

boundary 


=  0. 


(2.26) 


where  Xboundary  can  be  x,m„  or  xmax.  In  the  BTEC  1-D  Thermal  Model,  this  boundary  condition  is  implemented  as  an 
option  for  use  on  both  surface  boundaries. 


2.4.2  Insulating  Boundary  Condition 

The  insulating  boundary  condition  implies  that  no  energy  flows  through  a  boundary.  This  is  translated  mathematically 
to  imply  that  the  temperature  gradient  is  equal  to  zero  at  the  boundary: 


dv 

dx 


=  0 


Xboundary 


(2.27) 


2.4.3  Convective  Boundary  Condition 

The  convective  boundary  condition  is  a  linear  boundary  condition  for  which  the  rate  of  energy  loss  on  the  boundary  is 
proportional  to  the  temperature  difference  from  the  ambient  temperature.  The  proportionality  constant  has  a  value  in 
power  per  area  per  degree  and  is  known  as  the  convective  heat  transfer  rate,  he : 


dv 

dx 


=  he(v  -  Voo) 


dv 

dx 


=  ~he(v  -  Voo) 


Here  v  is  the  temperature  on  the  surface  and  vLX1  is  the  external  ambient  temperature. 


(2.28) 

(2.29) 


2.4.4  Radiative  Boundary  Condition 

A  radiative  boundary  condition  follows  a  Stephan-Boltzmann  law  where  the  energy  loss  rate  per  unit  area  is  propor¬ 
tional  to  the  difference  between  the  surface  and  ambient  temperatures  to  the  fourth  power: 

=  <re(T4  -  T4)  (2.30) 

Xmin 

T  is  in  units  of  Kelvin  and  Tr  is  the  effective  ambient  temperature.  The  effective  ambient  temperature  is  a  value  which 
is  slightly  larger  than  the  ambient  temperature  and  corrects  for  heated  air  from  the  surface  [5]. 

This  boundary  condition  is  a  non-linear  boundary  condition  which  cannot  easily  be  represented  in  a  finite  difference 
implementation.  The  following  Taylor  approximation  is  used  to  linearize  the  boundary  condition  [6]: 

T4  -  T4  ^  4T?(v  -  vr)  (2.31) 

Also,  the  heat  transfer  coefficient  can  be  defined  as 

hr  =  4ecrT  * .  (2.32) 
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The  linearized  radiation  boundary  condition  equation  can  be  written  as 


dv 

dx 


=  hr(v  -  vr) 


(2.33) 


and 


dv 

dx 


=  -hr(v  -  vr). 


(2.34) 


2.4.5  Evaporative  Boundary  Condition 


The  Lewis  analogy  for  evaporative  loss  is  implemented  in  the  model  as  an  optional  surface  boundary  condition.  The 
method  has  been  employed  to  include  the  energy  loss  from  surface  evaporation  as  a  function  of  ambient  temperature 
and  relative  humidity: 


dv 
C dx 

dv 
K  dx 


—  Qvap 
—  Qvap 


(2.35) 

(2.36) 


2.4.6  Combined  Boundary  Conditions 


As  these  surfaces  have  differing  signs  for  outward-facing  unit  normal  vectors,  energy  flows  will  be  described  with 
opposite  signs: 


dv 
C  dx 


he(v  -  Voo)  +  hr(v  -  vr)  +  Qv 


•ap 


dv 

(dx 


he(v  -  Voo)  +  hr(v  -  vr)  +  Q 


vap 


(2.37) 

(2.38) 


In  addition,  values  of  emissivity  and  convective  heat  transfer  rates  (lie  and  he)  may  have  differing  values,  if  material 
properties  differ  at  the  two  interfaces. 


2.4.7  Boundary  Condition  Derivations 

In  order  to  apply  the  boundary  conditions  to  the  Crank-Nicholson  Method,  the  derivatives  in  Equations  2.37  and  2.38 
need  to  be  approximated  with  a  finite  difference  and  then  solved  for  the  fictitious  points  v_i  and  vM.  This  needs  to  be 
done  for  vn  and  v"+1.  Lastly,  to  be  able  to  switch  any  of  the  three  boundary  conditions  on  or  off,  three  variables  a,  b, 
and  c  will  be  added: 


dv 

dx 


=  ahe(v  -  Vqo)  +  bhr(y  -  vr)  +  cQ 


vap 


Xmax 


To  find  vn_l  the  following  formulation  is  used: 


ahe(v  -  Voo)  +  bhr(y  -  vr)  +  cQmp 


(2.39) 

(2.40) 


kq 


X\  -  X_1 


=  ahe(vl  -  Voo)  +  bhri^Q  -  Vr)  +  c2,.flp(Vo) 


(2.41) 


Let 
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Pi  = 


Kj 

Ax, 


PI,  =  PI  -  —[ahe(pQ  -  Voo)  +  bhr{V'0  -  vr )  +  cQvap(VlQ)\ 
Pi) 


and  then  rearrange  terms  to  obtain 


1  1  c 

=  P\  +  —(alleVoo  +  bhrvr )  -  —  (ahe  +  bhr)vQ  -  —QvapiPo). 
P  0  P  0  P  0 


Finally  let 


qr  =  aheVco  +  bh,vr. 


qi  =  ahe  +  bhr. 


and 


v" ,  =  v?  +  ^  -  —  Gvb»(vS). 

-1  1  Po  Po  0  Po  0 


For  v^,  the  result  is  very  similar: 


To  find  v"^1  use 


P'm  -  P'm-2  +  n  n  P'm-  1  n  Qvapi^M-l ) 

Pa/-i  Pm-i  Pm-i 


y”+l  — 

*o~ - —  =  a/*e(Vo+1  -  Voo)  +  Wj,(vo+1  -  Vr)  +  cQvap{vl+l). 

X\  —  X-\ 

In  order  to  find  v^1  and  v^1,  the  value  of  QvapC*'”*1)  must  be  approximated  using  the  Taylor  series: 

n v  dQvclp 


Q,ap(P’+l )  *  GvaP(v") 


t/v 


(v"+1  -  v”) 


For  more  simplification,  let 


y  -  c- 


dQv, 


>cip 


dv 


By  putting  in  the  approximation  and  rearranging  terms. 


v"_;'  =  v^+1  -  —  [a/if(vg+1  -  Voo)  +  bhrivf1  -  vr)  +  cgVflj,(v£)  +  r(>o+I  -  vg)] 
Po 

can  be  found.  Next  use  qr  and  <72  to  once  more  simplify  the  answer: 

v"?  =  vr1  +  ^  +  £  -  -  -£)„(« 

■'  1  A)  A>  A>  A) 

Again,  the  v^1  point  is  very  similar: 

yi+i  _  yn+i  !  yvM-i  gr _ ®JlZv"+1 _ 0  fv"  ) 

yM  ~  yM-2  +  R  R  r  VM- 1  o  *!vaptyM-P 

Pm-i  Pm-i  Pm-i  Pm-i 


(2.42) 

(2.43) 

(2.44) 

(2.45) 

(2.46) 

(2.47) 

(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 

(2.53) 

(2.54) 
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2.5  Boundary  Condition  Implementation 

With  all  the  fictitious  end  points  found,  they  can  be  placed  into  the  Crank-Nicholson  Method. 

Start  with  equation  2.16  at  i  =  0: 

LHS,  =  flov"V  +  ^ovo+1  +  c^+1 
Incorporating  equation  2.53  into  equation  2.55  yields 

t  tic  _  ( tL  ^2  'y  L,\  n+l  .  (  L  .  ~L\  ti+\ 

LHS  o  —  u?o  —  a,  jv,  +  (c  q  +  a0)vi 

aL 

+  jQ{y<  +  qr  -  cQvap(v$). 

Note  that  the  last  term  in  equation  2.56  needs  to  be  moved  to  the  RHS  because  it  does  not  have  a  v"+1 
results  in  the  final  LHS  in  being 

LHS  o  =  | b(j  -  )vg+1  +  (Co  +  «oK+1- 

Start  with  equation  2.23  at  i  =  0: 

RHS o  -  +  bgvS  +  CqV[  +  2g, 

Putting  equation  2.47  into  equation  2.58  yields 

RHS,  =  {b*  -  g<)vS  +  (c*  +  al)V\ 

aR 

+  -  cQvap(V^)2g,. 

The  piece  from  the  LHS  is  added  on.  Then 


RHS,  =  (fc*  -  +  (a*  +  aR)V\ 

d L 

+  -j^{qr  -  cQvap(v,))  -  ^(r^o  +  q>-  -  cQvapW',))  +  ig. 


is  found.  For  M  -  1,  the  LHS  is 


LHS  M-l  -  (cM-l  +flM-l)VM-2  +  |^M-l  ^ 

and  the  RHS  is 

RHS  M-l  =  (CM- 1  +  aM-0VM-2  +  [b  M-l  ~  Pm  ~aM-ljVM-l 
+  JhrSqr  -  CQ''ap(V‘M-l))  ~  Pm  “  iy^M-  1  +  Qr  -  cQvap(y  M-  \ ))  +  ^ gO ■ 


(2.55) 

(2.56) 

term.  This 

(2.57) 

(2.58) 

(2.59) 

(2.60) 

(2.61) 

(2.62) 
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Chapter  3 

Two-Dimensional  Heat  Equation 


3.1  Heat  Equation 


Beginning  with  the  two-dimensional  heat  equation  in  cylindrical  coordinates. 


dv  Kdv  d  I  dv\ 
pC~dt  ~  ~r~d~r  +  Fr[Kd~rj 


+  A 


(3.1) 


where  k  =  k(z)  is  the  thermal  conductivity  of  the  tissue,  measured  in  [  CJS  c  ] ,  c  =  c(z)  is  the  specific  heat  of  the  tissue 
measured  in  p  =  p(z )  is  the  density  of  the  tissue  measured  in  |  |,  v  =  viz.,  r,  t)  is  the  temperature  rise  at  the 
coordinate  (z,  r,  t)  [°C],  and  A  =  A(z,  r,  t)  is  the  source  term  in  units  of  |  |. 


The  numerical  approach  used  in  the  BTEC  is  that  of  Peaceman  Rachford  [8]  and  utilizes  a  split  time  step  to  alter¬ 
nate  between  radial  and  axial  derivatives. 


3.2  Finite  Difference  Representation 

Simple  numerical  central -difference  derivative  approximations  are  shown  below  for  both  first  and  second  order  deriva¬ 
tives.  They  are  shown  first  for  a  uniform  grid  and  then  for  a  stretched  grid. 


3.2.1  Uniform  Grid 

First  Partial  Derivative 


Second  Partial  Derivative 


dv  _  vj+ 1  -  vh\ 
d  r  ~  2A  r 


d2v  _  vj+ 1  -  2  vj  +  Vj~  i 
dr2  ~  A  r2 


(3.2) 


(3.3) 


3.3  Non-Uniform  Grid 

In  order  to  extend  physical  space  out  in  the  radial  direction  without  using  vast  computational  space  or  making  the 
grid  spaces  about  the  axis  too  large  (where  the  highest  concentration  of  points  are  needed  in  order  to  model  a  beam 
accurately),  a  stretched  grid  is  required.  The  following  finite  differences  are  geared  to  handle  the  variable  step  sizes  in 
r.  The  notation  used  is  defined  below: 
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A r  =  rj+1  -  rhl 
A  r+  =  rj+ 1  -  rj 


3.3.1  First  Partial  Derivatives 


dy  =  o+ 1  ~  vj- 1 

<9r  A  r 


3.3.2  Second  Partial  Derivatives 


d2v 

dr2 


Q+  i  ~  VJ 
A  r+ 


yj-v7-i|  2 
Ar_  Ar 


(3.4) 


(3.5) 


3.3.3  Simplified  Notation 

The  finite  difference  approximation  for  the  numerical  solution  of  the  heat  equation  requires  the  following  terms  (from 
11.13):  v(z,  r,  t)  — »v?y- 

r— >r;j=0,l,2,...,N 

z  — >Zi  i=0,l,2,...,M 

A(z,r,t)  ->A'lj 

Where  i  is  the  axial  coordinate  (z)  index,  j  is  the  radial  coordinate  (r)  index,  and  n  is  the  time  index  representing 
the  time  step  At,,  from  an  initial  condition. 

The  derivation  will  be  simplified  further  by  using  the  following  defined  terms: 

P  =  Pl  =  2-w 


n  =  n  +  4 


n”  =  n  +  1 


A Ki  =  Ki+l  -  Ki- 1 

An  =  o+i  -  o-i 


A Zi  =  Zi+ 1  -  r,. 


A r+  =  O+i  -  O 


A  —  r,  —  r 


j  _  O-i 


Az+  =  Zi+1  -  Zi 


Az-  =  Zi  -  Zi 


i- 1 


=  /;Af„  (assuming  equally-spaced  time  steps) 
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3.3.4  Finite  Differences 

The  terms  required  to  represent  the  two-dimensional  heat  equation  in  a  finite  difference  form  become 


k  ldv\n 

r  W}u 

d  I  dv\n 

dr  \Kdrj  . . 


d_ 

dz 


n 

ij 


pc 


dv 

dt 


ij 


rArj 


^ ( vn  —  vn  y”  —  y”  \ 
2  Ki  viJ+ 1  Vi,j  Vi,j  Vi,j- 1 


Ar, 

2/t,- 

Azi 


Ar+ 


Ar_ 


tv", 

;+lj  ij  i,7  1-1,7 


Az+ 


Az_ 


A  «-</)• 


(3.6) 

(3.7) 


(3.8) 

(3.9) 


Axial  Boundary  Condition 

Due  to  the  symmetry  of  an  axial  problem 

dvn. 

-JT=0-  (3.10) 

dr 

By  evaluation  of  equation  3.10  using  the  finite  difference  representation  from  equation  3.6,  the  following  relationship 
can  be  found  across  the  r  —  0  boundary: 


<-t  =  Ki  (3-11) 

This  is  convenient  because  the  vni  x  point  doesn’t  exist  in  the  problem  space  and  is  just  considered  a  “phantom  point.” 
Also,  because  of  the  symmetry  of  the  problem,  the  phantom  grid  spaces  are  evaluated  by 


Ar_  =  Ar+ 
A  r  =  2A  r+  . 


(3.12) 


Applying  Equations  3.11  and  3.12  from  above  to  equation  3.7,  the  radial  second  derivative  finite  difference  approxi¬ 
mation  at  r  =  0  becomes 


d  I  dv\n 
d~r[KFr).  o 


2/e, 

A A 


yi,  l 


(3.13) 


The  relationship  in  3.10  does  not,  however,  imply  that  the  equation  3.6  goes  to  0  as  does  not  necessarily  equate  to  0. 
To  find  this  relationship,  L’Hospital’s  Rule  is  utilized: 


1  dv  d2' 


(3.14) 

Combining  equation  3.13  and  equation  3.14,  the  first  derivative  finite  difference  approximation  at  r  =  0  becomes 


/•->  0  r  dr  dr2 


K  tdvV1 

~r  Wi,o 


(3.15) 
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3.3.5  Finite  Differences  with  Simplified  Notation 

The  simplified  finite  difference  expressions  are  substituted  into  the  heat  Equation  3.1.  Note  that  there  will  be  a  separate 
expression  for  j  =  0  because  of  the  r  -  0  boundary  condition: 


j  =  0 


+ 


+ 


+ 


2k,  /y?+i,0-y?.o 

A z,  \  Az+ 


Az-  ) 


A" 

Ai,  0 


j  =  1.2,  .  ..N  -  1 


(3.16) 


fv".  ,  -  v"  ■  v".-v".  , '1 

/.j+i  i.j  _  hj  i,j- 1 

Ar+  Ar_ 

y»  _  y"  .  y"  .  _  y»  \ 

1+1 J  !J  1 J  1-1 J 


AO 

2k,- 

Az,- 


Az+ 


Az- 


—  (y«  .  _  y"  .) 

Vi+lj  Vi-Lj) 


A" 

AiJ 


(3.17) 


3.3.6  First  Half  of  Peaceman-Rachford 

The  Peaceman-Rachford  method  consists  of  a  split  time  step.  The  derivation  of  the  first  half  of  the  time  step  is  outlined 
in  this  section.  It  is  an  advantage  to  symbolically  represent  all  of  the  radial  differential  operators  as  D,  and  the  axial 
operators  as  Dz  so  the  problem  can  be  represented  in  a  familiar  matrix-vector  language. 

Equations  3.16  and  3.17  are  rearranged  to  yield 


j  =  0 
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tfXo-tfXo  - 


-i^Ly" 
A  ri  U  A ri  ifi 


2k, 


'  vn' 
y!+l,0 


Az,-Az+ 

5k  ,  6k  , 

+  A  zfi+1’°  A z2^1’0 

+  A"' 


2k, 


-VL,  n  - 


2a,- 


2a,' 

Az,-Az+  V'’°  +  Az,Az-  K,+1’°  AziAz-  V,’° 


-vfn 


2a,-  5a 

y"'  l  (  2Ki  l  2/C/  1 

J  + 

2a,-  5a  ' 

,  Az,-Az_  "  Az2, 

'  ''0  \Az,-Az+  AziAz-t 

*1,0  + 

,  Az,-Az+  Az2, 

v+i.o 


4a-,- 


4., 


5a:  2  a:,- 

Az2  AzjAz-  , 


<,,h-|«x.+|;(4  +  A-ih:0+ 


-W-£M£l*+A& 


5a  2a, 
k  Azr  Az,-Az+  , 


Ki+1,0 


(3.18) 


j  =  1,2,  ...N  -  1 


PKi-PKi  = 


— — y1!  . 


0A0  /'/A/;,  '•>  1 


2a,- 


-  AZi—S!  - 


2k‘  -yi+A^v" 


ArjAr+  iJ+1  Ar;Ar+  iJ  Ar;Ar_  iJ  ArjAr_  iJ  1 


2a,- 


V' 


2a,- 


-v? 


2a,- 


v?, 


2a,- 


-v?  , 


Az,-Az+  ,+1'2  Az,-Az+  I,j  Az,-Az_  !J  Az,-Az_  '  lj 
_5a  _  _5a 

+  Az2  i+1J  Az2  Mj‘ 


+ 


2a,-  5a 

/  ,  (  2k<  ,  ) 

v”'  4- 

2a,-  5a  ' 

,  Az,-Az_  Az2, 

'  Ij  \Az,-Az+  Az,-Az_ ) 

VU  + 

,  Az,- Az+  Az2 , 

#‘<1  + 


=  «V'+|-^-  +  ,  + 

Pl  «  l  TjATj  ArjAr_  1  ^ 


2a,- 


2a,- 


ArjAr+  ArjAr_  /  !J 


v?.-  + 


Ki+lJ 

2a,- 


ryAry  ArjAr+)  ,J+1  ,J 


+A", 


'  5a 

2a,-  ] 

'  5a 

2a,-  'j 

U? 

Az,-Az_ , 

^-ij  r  a^-  U-  az-//^ 

>1 

ft 

to  I 

Az,-Az+ , 

<+U 


/_AV_  2  A/ 

WAo  + 


<1-1 


—  (—  +  — ))y"+(  2*' 

Ary-  \  Ar+  Ar_ ))  ,J  \  ArjAr+ 


rAn 


v"  +  4" 

^/j+l  +  A/j 


(3.19) 
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In  order  to  simplify  the  matrix  representation,  the  following  terms  are  defined: 


Ko s  o 


drr , 


2  Kj  1 
A Zi  * 

[y  + 

y 

(A 

U> 

III 

^  6k 

{  y ' 

>  1  -£■ 
"•1  1 
+  to|~> 

) 

Oo 

J--  OJ 

o 

III 

I  4  Ki  \ 

W+) 

2  Kj  , 

[y  + 

1  V 

Oo 

OJ 

III 

(  2k, 

A  rj  1 

A  rj) 

^  A  rjAr+ 

-O 

rjArj  j 


Substitution  of  these  simplification  variables  into  3.18  and  3.19  yields 

v?-  i,o  +  1 v?+i,o  =  6rtoKo  +  +  A"o 


<t ,7  +  =  «j_,  +  SrlAj  +  +  4  ' 


In  matrix-vector  notation,  the  first  P.R.  step  becomes 

by'  =  DrV'  +  A"' 

com  —  by  +  a n 


6rlo 

K 

0 


OJifi 

OJi,  1 

W/,2 

. 

Wi,N-2 

.  Wi,N- 1  . 

6rlo 

0 

Sri,l 

dr3 

dr,1. 

dr3. 

6b 


uri,N-  2 


0 

Ar“ 

uri,N-2  uri,N-2 
uri,N- 1 


v/,o 

Ai.o 

v,-.i 

Ai,  l 

Vj,2 

Ai,2 

+ 

ViJV-2 

Ai,N- 2 

.  Aijf-i  _ 

Dzv"  =  w„ 


ft,  dz3  0 
&2  &2 

0  fej  fej  &3 


0 


& 


■2V— 2 

0 


6z 


M-2 


0 

&M-2 

^ZM-l 


y 

V2J 

U2J 

V2J 

— 

y.j 

VM-2J 

C>M-2,j 

-  VM-\,j  . 

.  C>M-l,j  . 

3.3.7  Second  Half  of  Peaceman-Rachford 

Now  equations  3.16  and  3.17  can  be  arranged  into  the  form  of  the  second  half  of  the  Peaceman-Rachford  step, 
done  in  the  same  manner  as  the  first  half  of  the  PR.  step. 


j  =  0 
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(3.20) 

(3.21) 


(3.22) 


This  is 


= 


± l y?"  -  ±Ly”" 
A r;  u  Ar;  ''-0 


1+1,0 


AziAz+ 

W* 

0  i+1,0 


Az? 


2*^.0 

Az,-Az+ 

*<-1,0 

Az? 


2^-v"  2k -v*1 

AK‘vi,  0  ^W-l.O 

+ 


AziAz- 


AziAz- 


A" 


i,0 


W  0  + 


^)v-;  +  ^v:  = 


f  +  / 
2k; 


Arr 


5k 


AziAz-  Azr 


-1.0 


2k,- 


2k, 


Az,-Az+  Az,'Az- 


Vf'n  + 


2k[ 

Az,-Az+ 


5k 

A? 


/>  +  4" 

4+1,0  +  Ai,0 


tf+^ko+(«= 


5k  2k, ■ 

Az;  Az/Az- 


AH 


1,2,  ...N  -  1 


Ibw1'  2k‘I  1 

„  1 1 

)<«* 

5k  2k,- 

r  az,  \az+ 

Az- , 

,Azr  '  AzAz+J 

v'1  4-  4'1 

4+i,o  +  A;,o 


(3.23) 


V’ 


0A0  «+1  r,Ar,  ■>  1 


2k, • 


« ,  - 


2k;  -v?,  _  2k, ■  ^ 


2k, • 


ArjAr+  U+1  Ary-Ar+  iJ  ArjAr.  u  Ar;Ar_  1 


2k, ■ 


-v?, 


2k,- 


2*'  >  . _ .  . _ .  .  T _  . 

Az,Az+  !+lj  Az,Az+  K,J  Az,Az_  V,J  ^  Az,Az_  lj 

5k  ,  5k  , 

-^vL,  -  — ?v, 


Az: 


,2  i+U 


Az 


2  4-1J 


py.i  + 


2k, ■ 


v” 


■I 


2k, ■  2k,' 

+ 


1  \ r; Ar;  ArjAr+)  Li  1  \ArJAr+  Ar;Ar_ /  !J  \  r7-Ar;-  ArjAr+j  !J+1 


2k, ■ 


v?,. 


2k, 

5k' 

v”'  4-  1 

2k, ■ 

2k, ■ 

\y  + 

2k, ■  5k  ' 

k  Az,-Az_ 

Az?J 

yi-ij  + 

Az,-Az+ 

AziAz- 

rj 

,Az;Az+  ’  Azr, 
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-  6%r)  <>-' + K +  %  (s- +  iz))  *&■ 

6k  2k,  \  ,  j  2k,  /I  1  \\  „>  6k  2ki  )  ,  , 

,  A zf  +  AziAz-J  ‘~hj  +  (  '  A zt  \Az+  +  Az.  ))  hj  +  Azj  +  Az„Az+ )  i+hj  +  ‘-J 

Again,  in  order  to  simplify  the  matrix  representation,  the  following  terms  are  defined: 

*J  =(-#  +  iSr)  &?  =  (e? -&(£  +  £))  &?  =(#  +  Efe) 


i,0  i, 0  y  i  Arl  I 

6r]  =(— - lKi  'l  6r2  =(/?"  +  —(—  +  — 

U  “  \riAri  ArjAr^j  ul  ij  ~  \Pi  T  A rj  \A r+  T  Ar_ 

Substitution  of  these  simplification  variables  into  3.23  and  3.24  yields 

Srlo^fi  +  Htffl  -  i.o  +  fc?y"o  + 


2*7  1 

A Zi  * 

[w+  + 

i)) 

6z]  =| 

(  6k  , 
[Azf  +  ‘ 

4Ki  ' 

) 

III 

O 
m  .«r 

°o 

[  4k‘  ) 

A  r\ 

I 

2k i  j 

f  1  + 

1  )) 

I 

f  lKi 

A  rj  ’ 

\ Ar+  + 

Arjj 

,  ArjAr. 

oAo  / 


+  ^u<o  =  <  i,o  +  &?v"o  +  *Mi,o  +  A"o 

+  +  *Xi,  +  <j  ■ 


In  matrix-vector  notation,  the  second  P.R.  step  becomes 


by"  =  by'  +  An' 


oj  =  by'  +  An' 


wi  j 
U2J 

60 

61 

&o 

6z] 

0 

0 

yW 

y2  j 

4«' 

Av 

A"  . 

2  j 

W3  j 

= 

0 

6z\ 

&2  *2 

0 

y3j 

+ 

A  n 

A3  J 

OJM-2,j 

0 

*^ZM-2  ^ZM-2 

0  <*4-t 

^ZM-2 

^4-1 

yM-2j 
-  yM-lj  . 

An' 

An' 

by  =  u> 


6rlo 

*y> 

0 

K 

5rh 

*4 

0 

5rll 

o 


6r 


l 

i,N-2 


0 


0 

Vi,0 

w;,  o 

yi,l 

w/,i 

Vj,:  2 

w;,2 

0 

dr* 

uri,N-  2 
Sr2 

uri,N- 1 

V,yv-1 

.  ViiAr_i 

.  t<4yv-i  . 

(3.24) 


(3.25) 

(3.26) 


(3.27) 
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3.4  Boundary  Conditions 

3.4.1  Constant-Temperature  Boundary  Condition 

Also  known  as  the  sink,  this  is  the  simplest  of  the  boundary  conditions.  This  boundary  condition  holds  a  boundary 
temperature  change  equal  to  zero.  In  the  model,  this  implies  the  following: 


v|w  =  0  (3.28) 

vL,  =  o 

y\Zmax  =  0  (3.29) 

In  the  implementation  of  the  BTEC  Thermal  Model,  this  boundary  condition  has  been  implemented  for  optional  use 
on  all  boundary  surfaces. 


3.4.2  Insulating  Boundary  Condition 

The  insulating  boundary  condition  implies  that  no  energy  flows  through  a  boundary.  This  is  represented  mathemati¬ 
cally  as  a  zero  temperature  gradient  at  the  boundary: 


dv 
dr  rn 
dv 
dz  Z: 
dv 


0 

0 

0 


Applying  a  central  finite  difference  to  equation  3.30,  it  is  found  that 


tf,N  ~  v! In-2 


where  rmax  corresponds  to  j=N-l. 


(3.30) 


(3.31) 


(3.32) 


The  first  half-time  step  in  the  finite  difference  Peaceman  Rachford  Method  (implicit  in  z)  is 


dz, 


1 ,  n' 


■  &?v" 


-i  +  ^i^+ijv-i  _  driN_ iViN_2  ■ 


dr, 


iJV-l  vi,N- 1 


’  drhN_  |  v"N  +  AiJS[_ j  . 


(3.33) 


Incorporating  equation  3.32  into  3.33,  the  following  relationship  is  derived  to  be  implemented  for  an  insulating  bound¬ 
ary  condition: 


i,n-i  +  s^n_ j  +  =  (drjj  +  Sijj)  v"iN_2  +  6 +  A* 


(3.34) 


The  second  half-time  step  of  the  Peaceman  Rachford  Method  (implicit  in  r)  has  the  finite  difference  representation 
shown  below: 


dr 


1  vn" 
i,N-l  vi,N- 2 


dr 


i.N-1  vi,N- 1 


•  i.jv-1  +  dzfv&r-i  + 


+  A " 

i+1  ,iV —  1  ^  Ai,N- 


(3.35) 


Applying  equation  3.32  into  3.35,  the  following  relationship  is  derived  to  be  implemented  for  an  insulating  boundary 
condition: 


(^ri'.JV-l  +  )rf,N-2  +  drfjt-irfjf-i  ~  dz]V}_  1  A 


-l  +  dZj  v?N_  i 


+  An 

-1  1 


(3.36) 
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3.4.3  Convective  Boundary  Condition 


The  convective  boundary  condition  is  a  linear  boundary  condition  for  which  the  rate  of  energy  loss  on  the  boundary 
is  proportional  to  the  temperature  difference  from  the  ambient  temperature.  The  proportionality  constant  has  a  value 
in  power  per  area  per  degree  and  is  known  as  the  convective  heat  transfer  rate,  he.  In  the  case  of  the  2-D  model,  this 
boundary  condition  has  been  implemented  as  an  optional  boundary  condition  for  the  z.min  and  z„,ax  surfaces: 


dv 

% 


=  he(v  -  Voo) 


Here,  v™  is  the  ambient  temperature,  relative  to  the  baseline  tissue  temperature. 


(3.37) 


3.4.4  Radiative  Boundary  Condition 

A  radiative  boundary  condition  follows  the  Stephan-Boltzmann  law  where  the  energy  loss  rate  per  unit  area  is  propor¬ 
tional  to  the  difference  between  the  absolute  surface  and  ambient  temperatures  to  the  fourth  power: 

=  cre(04  -  0?)  (3.38) 

The  temperatures,  6  and  9r,  are  the  surface  temperature  and  effective  ambient  temperature,  in  Kelvin.  The  effective 
ambient  temperature  is  slightly  higher  than  the  actual  ambient  temperature,  and  corrects  for  heated  air  from  the  surface. 
This  boundary  condition  is  a  non-linear  boundary  condition  which  cannot  easily  be  represented  in  a  finite  difference 
implementation. 


3.4.5  Evaporative  Boundary  Condition 


The  Lewis  analogy  for  evaporative  loss  is  implemented  in  the  model  as  an  optional  surface  boundary  condition.  The 
method  has  been  employed  to  include  the  energy  loss  from  surface  evaporation  as  a  function  of  ambient  temperature 
and  relative  humidity: 


=  Q 


vap 


(3.39) 


3.4.6  Combined  Boundary  Conditions 

In  order  to  implement  a  meaningful  surface  boundary  condition  for  the  2-D  model  in  cylindrical  coordinates,  it  is 
important  to  place  surfaces  at  the  minimum  and  maximum  z  coordinates.  As  these  surfaces  have  differing  signs  for 
outward-facing  unit  normal  vectors,  energy  flows  will  be  described  with  opposite  signs: 

=  he(v  -  vra)  +  creidf4  -  0?)  +  Qvap  (3.40) 


=  he(v  -  Voo)  +  ce(04  -  0?)  +  Qmp  (3.41) 

In  addition,  values  of  emissivity  and  convective  heat  transfer  rates  may  have  differing  values  if  material  properties 
differ  at  the  two  interfaces.  The  implementation  at  the  zmin  surface  within  the  finite  difference  method  are  considered 
and  then  the  additional  parameters  required  for  code  implementation  are  examined. 
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3.5  Boundary  Condition  Implementation 


In  the  BTEC  (2-D)  Peaceman  Rachford  Method,  there  is  an  alternating-direction  implicit  method.  For  each  half  time 
step,  the  system  of  equations  implementing  boundary  conditions  is  different.  The  half-step  implicit  in  z  is  considered 
first.  The  following  equation  utilizes  the  notation  for  differential  operators. 


SzjvUj  +  SzXj 


■  =  5rhvlj- 1  +  6rf,jVL 


(3.42) 


The  Zmin  coordinate  is  indexed  in  the  model  as  i  =  0,  such  that  the  surface  boundary  condition  can  be  represented  as 


*0^1, 


+  6z20vn0j  +  dzlv'lj  =  j  +  5rl/0j  +  5rl/Q 


+  A" 

1+1  +  A0J  ■ 


(3.43) 


The  v"  j  .  coordinate  is  a  fictitious  point,  eliminated  by  solving  for  this  coordinate  in  the  boundary  condition  represen¬ 
tation: 


/  /I  ,  dr  Qvap 

y- U  -  v'lj  +  p(hl+  hr)vn0J  +J  + - p 


(y'o  j)  r  (vg..)  , 


p 


(3.44) 


The  fictitious  point  is  eliminated  through  substitution  to  eliminate  vn_[  .  in  the  finite  difference  representation.  The  LHS 
of  the  finite  difference  representation  at  the  boundary  becomes 


LHS  -  vSj 


I  1  r(vo.P 

8z0(-(hi  +  hr)  +  — — — 


)  +  &, 


+&o 


qr  Qvap(y  Qj) 


p 


p 


+  V(j 
rO'oj) 


■&0 


p 


(3.45) 


The  second  half-time  step  of  the  Peaceman  Rachford  Method  (implicit  in  r)  has  the  finite  difference  representation 
shown  below: 


K/C  ,  +  Sr^  +  Prfjrfj+i  =  +  SzfV'j  +  +  A"' 

At  the  minimum  z  coordinate  boundary,  i  =  0,  yielding 

Srojvoj-i  +  6roA,j  +  6rojvoj+i  ~  &ov-ij  +  ^ovo,;  +  &0ylj  +^0j' 
and  the  explicit  boundary  condition  representation  is 


(3.46) 


(3.47) 


v-ij  -  v",j  +  +  brKj  +  ^  +  pQ'vpW 3j)-  (3.48) 

The  boundary  condition  relationship  in  equation  3.48  is  used  to  replace  the  phantom  point  in  the  finite  difference 
representations.  The  right  hand  side  of  second  PR.  step  then  becomes 


RHS  =  vgj. 


Szl0-(hi  +  hr)  +  Szl 


[&i  +  &o] 


+&n 


^  +  pQvap(VQj) 


+  A ” 

+  A0 J  ■ 


(3.49) 


The  implementation  of  surface  boundary  conditions  involve  several  functions  determined  by  absolute  temperature  in 
Kelvin,  temperature  in  Celsius,  or  a  relative  temperature.  The  finite  difference  solution  is  implemented  in  terms  of 
change  in  temperature.  The  code  tracks  various  temperatures  in  terms  of  the  current  finite  difference  solution  or  in 
terms  of  measured  temperature. 
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In  addition,  the  software  implementation  of  the  surface  boundary  conditions  will  allow  the  user  to  simply  switch  a 
given  effect  on  or  off.  The  code  will  use  the  same  solution  for  any  combination  of  convective,  radiative,  and  evapora¬ 
tive  effects. 

Various  temperatures  are  defined  to  simplify  the  notation  of  the  boundary  conditions.  Temperatures  in  Kelvin  are 
denoted  by,  6 ,  temperatures  in  Celsius  are  denoted  by,  T .  The  convention  of  relative  temperature  has  already  been 
used  throughout  this  paper,  and  are  denoted  by  v. 


v=T -T0 


(3.50) 


where  To  is  the  baseline  tissue  temperature.  The  relative  ambient  temperature  is 


Voo  =  Tm  -  To  . 


(3.51) 


The  effective  relative  ambient  temperature  is 

v,.  =  Tr  -  To  (3.52) 

In  addition,  the  variables  a ,  b,  and  c  are  debited  to  represent  switch  values  for  each  of  the  convective,  radiative,  and 
evaporative  boundary  conditions,  respectively.  The  following  intermediate  values  are  also  used  to  simplify  notation. 


hr  =  AecrT ,3 

(3.53) 

qi  =  ah  i  +  bhr 

(3.54) 

=  ah  iVoo  +  bhrvr 

(3.55) 

dQvap  | 

'  dT  'v 

(3.56) 

The  boundary  condition  at  the  zmin  surface  is  then  given  by 


=  ahr(v  -  Voo)  +  bcre(ff*  -  (fr)  +  cQvap  ,  (3.57) 

where  a ,  b,  and  c  may  be  assigned  differing  values  on  the  zmax  surface. 

The  derivation  of  the  implicit  and  explicit  (in  z)  bnite  difference  representation  of  the  i  =  0,  or  zmin  grid  positions  are 
obtained  in  the  following  equations : 


v? 


0  J 


c  2  -i,92  y 

&o  ~  6zq(j  +  p) 


+  </(&0  +  &o) 


+&n 


-  fcojVgjL.  1  +  SfqjVqj  +  +  Aqj 


Sro./a,j-\  +  6ro,/oj  +  <^o,/oj+i 


K  -  6zo(-p  +  -p) 


+  V(j(6z'o  +  &o) 


+Sz’ 


q-T--0  (vn  ) 

P  UvapKYoj) 


+  A " 

+  A0,j 


(3.58) 


(3.59) 
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The  same  boundary  conditions  at  the  zmax  surface  are  also  defined. 

Note:  In  the  BTEC  model  implementation,  the  nz  -  1  index  is  a  position  on  the  boundary  and  the  nz  index  repre¬ 
sents  a  fictional  or  phantom  grid-point,  and  repeats  the  derivation  for  the  equations  above: 


v'mj 


5zm  -  SzU 


~3l  +  Zl 
P  P 


+  VM-lJ^ZM  +  &M) 


+8z 


M 


p  pQ^p^hj)  + 


-  SrM,jVM,j-\  +  ^rM,jVM,j  +  drM,jVM,j+ 1  +  AM,j 


^rM,jVM,j-i  +  +  drM,jVM,j+ 1 


- 


j) 


+  Pm-ij^Zm  +  8z3M) 


+8zl 


^--0  (v"  ) 

P  p'*lvap\v  M,j> 


4-  A" 

+  AM,j 


(3.60) 


(3.61) 
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Chapter  4 


Source  Terms  and  Perfusion  Effects 


4.1  Source  Term  Representation 


The  BTEC  thermal  model  has  the  capability  of  using  several  different  types  of  source.  Ultimately,  a  3-dimensional, 
cylindrically  symmetric,  time  dependent,  source  term  with  dimensions  of  W/cm 3  is  needed  to  solve  the  heat  equation. 
All  source  terms  (from  multiple  lasers,  perfusion,  etc.)  can  be  summed  into  one  function,  denoted  by  A (z,  r,  t).  This 
function  is  then  the  single  source  term  for  the  Heat  Equation. 

BTEC  thermal  can  compute  a  Beer’s  Law  source  based  on  the  profile  of  the  beam  and  the  absorption.  With  the  Beer’s 
Law  source  the  model  is  set  up  to  compute  the  source  term  for  a  flat-top,  Gaussian,  and  annular  profile  beams.  The 
BTEC  model  can  also  utilize  wave  propagation.  The  finite  difference  wave  propagation  and  Hankel  wave  propagation 
methods  compute  the  source  term  based  on  the  initial  electric  field  distribution  and  the  current  refractive  index  of  each 
node.  The  refractive  index  varies  with  temperature  and  the  source  term  is  recomputed.  The  wave  propagation  methods 
and  their  derivations  are  discussed  in  more  detail  in  Chapter  6. 

The  heat  source  term  for  the  formulation  of  the  thermal  model  was  examined  through  two  differing  approaches.  The 
different  approaches  are  referred  to  as  the  optical  model  and  the  radio  frequency  model.  In  the  optical  effects  formu¬ 
lation,  the  time-dependent  solution  to  the  Heat  Equation  is  determined  for  a  source  term  debited  by  equation  4.5.  A 
localized  heat  source  term  estimated  from  a  simple  linear  absorption  is  commonly  used  for  optical  frequencies.  The 
time  dependent  amplitude  is  assumed  to  be  a  square-pulse  temporal  prohle  with  a  Gaussian,  top-hat,  or  annular  spatial 
distribution. 


4.1.1  Beer’s  Law 


The  3-dimensional,  time  dependent,  source  term  computed  by  Beer’s  Law  is 


A(z,  r,  t)  =  /uMr,  t)e  ^ 


(4.1) 


The  exponential  decay  in  the  z  direction  is  dependent  on  the  absorption  coefficient,  \xa.  7o(r,  t)  is  the  incident,  spacial 
beam  prohle.  The  three  beam  profiles  mentioned  are  represented  below. 

Flat-Top 


(4.2) 


Gaussian 


Io(r,t)  = 


(4.3) 
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Annular 


h (r,  t )  = 


Iq,  CT\  <  r  <(T2 
0,  r  <  ct\,  r  >  <X2 


(4.4) 


The  beam  radius  is  given  by  , sigma,  and  the  Peak  Incident  Irradiance  is  Iq. 

4.1.2  SAR 


A(z,  r,  t )  =  h(z,  r)H0(A,  t)/ja(z.  A)  (4.5) 

This  source  term  provides  a  time-dependent  description  of  the  linear  absorption  of  optical  energy  as  a  function  of 
depth  in  the  tissue,  complete  with  spectral  and  radial  dependence  of  energy  being  absorbed.  The  variable  A  refers 
to  the  wavelength  of  the  THz  source.  The  function  h(z,  r)  specifies  the  relative  irradiance  for  a  given  position  in  the 
cylindrical  coordinate  system  and  includes  losses  due  to  linear  absorption.  This  term  also  addresses  the  focusing  of 
the  beam  through  the  tissues  using  a  specified  beam  waist  location  and  a  hyperbolic  function  to  assign  beam  radius 
as  a  function  of  position.  The  function  Hq(A,  t )  provides  the  maximum  irradiance  per  wavelength  division  at  a  given 
time  [W/crrr/nm],  The  value  of  pa{z,A)  represents  the  absorption  coefficient  [1/cm]  at  a  given  wavelength  within  the 
tissue,  which  is  determined  by  the  tissue  type  at  the  given  axial  depth,  z  [cm]. 

In  the  radio  frequency  formulation  of  the  source  term,  a  one-dimensional  Finite  Difference  Time  Domain  (FDTD) 
method  is  used  to  predict  the  specific  absorption  rate  (SAR)  within  single-  and  multi-layer  homogeneous  skin  slab 
models.  The  FDTD  method  is  an  explicit  time-domain  numerical  approach  for  solving  Maxwell’s  equations  in  a  dis¬ 
cretized  space[14,  11].  This  method  has  become  widely  used  for  predicting  the  electromagnetic  fields  and  energy 
absorption  rates  within  digital  voxelized  models.  Equations  4.6-4. 7  below  give  the  FDTD  update  equations  for  the  Ex 
and  Hy  fields,  respectively: 


i  C<  Af 

1  2ek 

EAnr-+ 

At 

ek 

1  +  9* 

2ek 

1  +  2* 

2  €k 

Az 

H. 


in+l_ 


yi  k 


1  <r‘Af 

1  2/4 

A t_ 

Hk 

1  +  -J— 

2/4 

1  + 

2/4 

M 


Ll<t+A 


Az 


(4.6) 


(4.7) 


In  Equations  4.6-4. 7,  cr  is  the  conductivity  in  Siemens  per  meter,  p  is  the  permeability  in  Henrys  per  meter,  and  cr* 
is  the  magnetic  loss  term  in  ohms  per  meter.  The  magnetic  loss  is  assumed  to  be  zero  and  the  localized  SAR  values 
within  the  tissue  are  calculated  by  equation  4.8.  The  index  k  represents  the  discretized  axial  coordinate,  z  =  k(Az),  and 
n  indicates  the  time  step,  t  =  n(At).  The  SAR  value  is  subsequently  defined  from  the  electric  field  strength  and  the 
density  of  the  tissue,  as  given  by  equation  4.8: 


SAR  = 


cr\E\2 

P 


In  equation  4.8,  p  is  the  density  [kg/cm3],  and  E  is  the  root-mean-square  of  the  electric  field. 


(4.8) 


4.2  Perfusion  Effects 

The  BTEC  thermal  model  employs  a  non-localized,  uniform  perfusion  term  to  simulate  the  convective  loss  of  energy 
due  to  blood  flow.  The  term  physically  represents  a  uniform  constant  exchange  of  mass  at  each  tissue  location,  with  the 
entering  mass  having  a  temperature  of  the  reference  tissue  temperature  (body  temperature).  The  form  of  the  perfusion 
source  term  equation  is  given  below: 


q(z,  r,  t )  =  — wp(z )  c(z )  v(z,  r,  t )  (4.9) 
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The  perfusion  flow-rate  term,  wp,  expressed  in  W/cm3 /s,  is  assigned  to  each  tissue  layer  type  within  the  model  and 
added  to  the  current  source  term  calculated  from  the  emitter’s  properties.  The  values  c  and  v  represent  specific  heat 
and  temperature  rise  in  the  tissue,  respectively.  The  perfusion  term,  c(z),  currently  implemented  is  constant  with 
temperature.  Typical  values  can  be  found  in  the  text  by  Welch[13], 
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Chapter  5 

Gaussian  Beam  Propagation 


Two  forms  of  propagation  analysis  are  currently  employed  in  the  BTEC  model:  a  Gaussian  Beam  Propagation  Method 
(GBP)  and  wave  propagation.  The  GBP  allows  for  Gaussian  beams  to  be  translated  through  linear  optics.  It  is  mostly 
used  to  compute  the  initial  conditions  for  the  wave  propagators  employed  in  the  BTEC  project.  For  the  Gaussian  Beam 
Propagator  in  this  project,  a  general  expression  using,  the  ABCD  matrix  method  common  in  optical  transformations  is 
needed  to  calculate  how  a  laser  beam  propagates  through  linear  geometric  optics.  Many  books  and  articles  solve  this 
problem  with  the  condition  (or  assumption)  that  the  lens  or  optical  component  is  placed  at  the  waist  of  the  approaching 
beam.  The  goal  is  an  expression  for  an  optical  system  placed  in  any  arbitrary  position  along  the  approaching  beam  path. 

The  electric  field  E  for  a  Gaussian  beam  can  be  defined  as 


C  A  w° 

t  =  A — exp 
0) 


+  l 


k(z  -  Za> )  -  tan 


(Z-Zu] 

kr2 

kr4  1) 

l  *  J 

2R(z)_ 

[8R3(z)Jj 

(5.1) 


where  a>o  is  the  beam  diameter  at  the  beam  waist  and  oj(z)  is  the  beam  diameter  at  a  z[cm]  coordinate  defined  by 


u(z)  =  O)0 


(5.2) 


O)  =  0)0 


(5.3) 


The  Raleigh  range  zo  can  be  define  as 


(5.4) 


(5.5) 


The  BTEC  model  follows  the  Milonni  and  Eberlys  LASERS  [4].  The  “Entry”  matrix  in  figure  5.2  is  added  to  account 
for  the  distance  from  the  beam  waist  that  the  geometric  optic  is  placed.  This  allows  for  optics  to  be  placed  at  any 
arbitrary  position  along  the  optical  axis.  This  figure  also  shows  matrix  multiplication  is  evaluated  from  right  to  left. 


29 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


A 

C 


Figure  5.1:  Matrix  Thin  Lens  Formulation 

Exit  Thin  Lens 


Entry 


1 

d' 

1 

0 

1 

d 

0 

1 

X 

-1/f 

1 

X 

0 

1 

— 

\” 

\ 

— 

y 

■l 

1 

d' 

1 

d 

N 

0 

1 

X 

-1/f 

-d/f+1 

A 

B 

_ 

l-d'/f 

d+d' 

C 

D 

-1/f 

1 

Figure  5.2:  ABCD  Matrix  Breakdown 


The  components  for  ABCD  matrix  are  placed  into  the  equation 

Aqt  +  B 


clf  = 


(5.6) 


-Cqi  +  D 

and  then  it  is  set  equal  to  the  definition  of  qj.  Equation  5.6,  and  the  real  and  imaginary  parts  of  both  sides  are  equated: 

1  1  (5.7) 


Qf  i  ,  _u_  i  .  a 

R  TUD2  R  TCCO2 


1  12  1 


q(z)  R(z )  ik  u>2(z) 

k =* 


(5.8) 


(5.9) 
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The  new  beam  waist  can  be  found  by  determining  where  the  radius  of  curvature  of  the  beam  approaches  infinity.  By 
substituting  this  value  back  into  the  equation  for  w(d),  the  relationship  for  a>o  can  be  found: 


a>0  -  COq 


_ Zr _ 

(Cz0-O)2+(C24) 


ZR 


Once  the  new  spot  size  is  known,  the  new  Rayleigh  range  can  easily  be  obtained: 


(5.10) 


-(Azo  -  B)(Cz0  -D)  +  ( ACzr ) 
0 Czo  -  D)2  +  (Ch2R) 


(5.11) 


^-^4)  ,5-12) 

By  finding  a  general  expression  for  the  new  beam  waist  and  Rayleigh  range,  different  optical  elements  can  be  accounted 
for  by  using  their  ABCD  values: 
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Figure  5.3:  Common  ABCD  Matrix  Representations 

Several  examples  of  common  usages  of  this  ABCD  method  are  shown  in  figure  5.3  and  their  parameters  are  outlined 
below: 

Planar  Bound  Propagation 


A  =  1 
B  =  0 
C  =  0 

)  =  ^i 
«2 


Thin  Lens  Propagation 


A  =  1 
B  -  0 


D  =  1 


Spherical  Lens  Propagation 
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A  =  1 
B  -  0 

112  ~  Hi 

n2Rc 
D  =  1 


where  R,  is  the  radius  of  curvature  of  the  lens. 
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Chapter  6 


Two-Dimensional  Helmholtz  Wave 
Equation 

6.1  Derivation  of  the  Theory 

It  is  important,before  looking  into  the  modeling  methods  that  are  used,  to  review  how  the  theoretical  equations  were 
found  and  what  approximations  were  used.  This  chapter  includes  a  derivation  of  the  equation  used  in  the  FD-BPM 
(Finite  Difference  Beam  Propagation  Model)  starting  from  Maxwell’s  equations. 


6.1.1  Maxwell  Equations 

The  derivation  of  the  wave  propagation  equation  used  in  the  BTEC  model  starts  with  the  differential  forms  of  Maxwell’s 
equations : 


V  •  E  =  0 

(6.1) 

<1 

w 

II 

o 

(6.2) 

<9H 

VxE^- 

(6.3) 

<9E 

V  x  H  =  e— 

(6.4) 

There  is  assumed  to  be  no  free  current.  Permittivity  and  permeability  are  assumed  to  be  independent  of  time.  This 
is  true  on  the  time  scale  of  light  propagation.  However,  these  properties  can  change  on  the  time  scale  of  the  thermal 
model  and  are  taken  into  account  in  the  changing  refractive  index. 

To  derive  the  Wave  Equation,  the  curl  of  Faraday’s  equation  (6.3)  is  taken: 

VxVxE  =  -Vx  yU—  (6.5) 

ot 

V  x  V  x  A  =  V(V  ■  A)  -  V2A  (6.6) 

Using  the  vector  identity  6.6  and  equation  6.5: 
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A  C\ 

V(V- ET-  V2E  =  -u-HV  x  H)  (6.7) 

at 

The  divergence  of  E  goes  to  zero  from  equation  6.1.  Substituting  equation  6.4  for  the  curl  of  H,  the  Wave  Equation  is 
found. 


6.1.2  Wave  Equation 


V2E(r,  t) 


<92E(r,  f) 


This  is  the  Wave  Equation.  From  separation  of  variables  a  solution  of  the  form 


E(r,  t)  =  G(r)e^' 

can  be  assumed. 


(6.8) 


(6.9) 


6.1.3  Helmholtz  Equation 

By  substituting  equation  6.9  into  equation  6.8, 

V2£(r)  +  yu«u2£(r)  =  0  (6.10) 

is  found.  This  is  the  Vector  Helmholtz  Equation  in  an  arbitrary  spatial  coordinate  system.  It  is  no  longer  time  depen¬ 
dent. 


The  Scalar  Helmholtz  Equation  is 


V2£(r)  +  yuew2£(r)  =  0  . 


In  cylindrical  coordinates  the  Laplacian  is 


0 


1  d  I  dA 


V  A  =  \r—  +  - 


1  dy 


d2A 


(6.11) 


(6.12) 


r  dr  \  dr  J  rjf  deft1  dz 2 

Assuming  a  problem  space  with  azimuthal  symmetry,  the  0  derivative  can  be  ignored. 

Also,  =  k2  —  k^n2  where  k(]  is  the  free  space  wave  number  and  n  is  the  refractive  index  of  the  mate¬ 

rial.  Using  this  and  equation  6.12,  the  Helmholtz  equation  6.11  in  cylindrical  coordinates  is 


d2&  1  (dS 

dz 2  r 


dr 


d2S  o  , 

+  +  V2£  =  o. 


A  plane  wave  solution  of  the  form 

6(r,z )  =  T'(r,z)<T,w 

is  now  assumed,  where  na  is  the  reference  index. 


(6.13) 


(6.14) 


Inserting  the  assumed  solution  6.14  into  the  Helmholtz  equation  6.13,  the  following  equation  is  derived: 


d2xV  d'V 


d2xV  1  d'V 
dr2  r  dr 


+  k20[n2(r,z)-n20\'¥ 


(6.15) 


equation  6.15  is  used  in  two  different  solutions.  The  first  is  a  paraxial  approximation  which  assumes  there  is  no  strong 
focusing  or  defocusing  effects.  It  utilizes  an  approximation  of  an  exponential  operator  as  a  rational  operator.  This 
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approximation  is  known  as  the  Caley  approximation.  [9]  The  second  numerical  solution  allows  for  strong  focusing 
effects  and  utilizes  Pade  approximate  operators.  Both  are  based  off  of  equation  6.15  and  will  be  discussed  in  greater 
detail  later. 


6.2  Numerical  Approach 

6.2.1  Finite  Differences 

The  following  finite  difference  approximations  are  used  in  the  numerical  derivations.  Finite  differences  allow  the 
representation  of  continuous  differential  equations  in  discrete  computational  space.  The  finite-  differences  are  shown 
below  for  both  first  and  second  order  derivatives.  They  are  shown  first  for  a  uniform  grid  and  then  for  a  stretched  grid. 


Uniform  Grid 

First  Order 


Second  Order 


a*  _  'Vj+i-'Vj-i 

dr  2A  r 

JV  =  T..,  -2‘1';~T,  , 

dr2  A  r2 


(6.16) 

(6.17) 


Non-Uniform  (Stretched)  Grid 

In  order  to  extend  the  physical  space  out  into  the  radial  direction  without  using  vast  computational  space  or  making 
the  grid  spaces  about  the  axis  too  large  (where  the  highest  concentration  of  points  are  needed  in  order  to  model  a  beam 
accurately),  a  stretched  grid  is  required  [3].  The  following  finite  differences  are  geared  to  handle  the  variable  spacial 
step  sizes.  The  notation  used  is  defined  below: 


First  Order 


Second  Order 


A r  =  rj+ 1  -  r>i 
A  r+  =  r j+]  -  rj 

Ar-  =  n  -  n- 1 


i 

+ 

I 

r 

(6.18) 

dr 

Ar 

d2xv 

Vj  ~  ^./-i  1 

2 

(6.19) 

dr 2 

Ar+ 

Ar_ 

Ar 

6.2.2  Boundary  Conditions 

These  boundary  conditions  apply  to  both  the  Caley  method  and  the  Pade  method.  These  methods  exploit  the  symmetry 
of  a  laser  beam  to  compute  a  solution  at  the  r  =  0  boundary.  They  also  allow  circumvention  of  the  division  by  0  in 
equation  6.15. 

Due  to  the  symmetry  of  an  axial  laser  beam. 


ff¥(r  =  0,z) 
Jr 


=  0  . 


(6.20) 


By  evaluation  of  equation  6.20  using  the  finite  difference  representation  from  equation  6.18,  the  following  relation  can 
be  found  across  the  r  —  0  boundary: 
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(6.21) 


f  ,  =  VF', 


This  is  convenient  because  the  point  doesn’t  exist  in  the  problem  space  and  is  considered  a  “phantom  point.” 
The  relation  in  6.20  does  not,  however,  imply  that  |  =  0  as  [j  does  not  necessarily  equate  to  0.  To  find  this 

relationship,  L’Hospital’s  Rule  is  utilized. 


1  ffV  _  d2xV 

r-> o  r  dr  dr 2 

Also,  because  of  the  symmetry  of  the  problem,  the  phantom  grid  spaces  are  evaluated  by 


(6.22) 


Ar_  =  Ar+ 
A  r  =  2Ar+ 


(6.23) 


6.3  Caley  Method 

The  method  of  interest  for  use  in  the  BTEC  thermal  model  is  the  Pade  Method.  The  Caley  Method  was  implemented 
first  because  it  was  conceptually  easier  to  understand.  However,  the  Pade  Method  is  a  more  powerful  method  as  it 
lacks  the  paraxial  assumption.  The  implementation  of  the  Caley  Method  gave  results  against  which  the  Pade  Method 
implementation  could  be  validated. 

In  the  Caley  form  solution,  a  paraxial  approximation  is  made  to  the  wave  equation.  If  it  is  assumed  that  there  are 
no  strong  focusing  or  defocusing  effects  occurring  in  the  media,  the  second  axial  derivative  can  be  assumed  to  be 
negligible  and  dropped  [2]: 


-likono 


d'V 

dz 


d2xV  1  d'V 
dr2  r  dr 


kl  [n2(r,z)  -  Hq]  T 


(6.24) 


The  equation  is  rearranged  for  the  ease  of  calculation  to  yield 


.d'V 

‘~dz 


-1 

2kriQ 


&_  \d_ 

dr2  r  dr 


\n2(r,z)  ~  «o| 


'F 


(6.25) 


i-r-  =  SV(z)  •  (6.26) 

dz 

In  equation  6.26,  an  operator  has  been  defined  for  ease  of  calculation.  This  S  will  contain  the  derivatives  in  the  r 
direction  as  well  as  constants.  It  is  defined  as 


S 


-1 

2knQ 


cP_  ]_d_ 
dr2  r  dr 


kl  [n2(r,z)  -  «q] 


(6.27) 


The  solution  to  equation  6.26  is  found  by  separating  and  integrating.  With  appropriately  chosen  bounds,  the  solution 
appears  as 


*F(r,z  +  Az)  =  e-lAzS'V(r,z)  . 


(6.28) 
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6.3.1  Caley’s  Form  of  the  Exponential  Operator 


The  Caley  method  as  it  is  referred  to  here  is  a  Crank-Nicholson  Method  that  uses  the  Caley  form  of  the  exponential 
operator.  The  Caley  form  is  an  approximation  that  allows  the  exponential  operator  to  be  represented  as  a  rational 
operator: 


e~iSAz 


1  -  }iSAZ 
1  +  {iSAz 


Using  the  Caley  relationship  6.29,  equation  6.28  becomes 


(6.29) 


,  iS  A  z 

,  iS  Az 

1  + 

2 

+(z  +  Az)  = 

1 

2 

'Viz). 


(6.30) 


Substituting  the  transverse  derivatives  and  constants  contained  in  the  S ,  from  equation  6.27  and  the  finite  difference 
representations  of  these  transverse  derivatives  from  Equations  6.18  and  6.19,  equation  6.30  becomes  Equations  6.31 
and  6.32.  Equation  6.31  is  the  left  side  of  equation  6.30  and  6.32  is  the  right  side: 


LHS 
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1  4k0n0 


/VI/I+1  _  11/1+ 1  VT//+1  _  U//+1  \ 

7+1  J  J  7-1 


Ar+ 


A  r- 


-  vpi+1  _  Vf/i+1 

2  ,  1  7+1  7-1  ,  ,2  f  2 


A  r  r  A  r 


kl[n2(r,z)  -  nl]  +' 


(6.31) 


RHS 
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(6.32) 


where  the  following  notation  is  again  adapted: 


T(r,z)  =  T' 


*V(r  -  A r,z  +  A z)  =  'F'+11 


6.3.2  System  of  Equations 

In  order  to  easily  view  the  system  of  equations  formed  in  6.31  and  6.32,  the  equations  are  ordered  by  the  T*  values  at 
each  point  in  the  r  and  z  directions. 

LHS 
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(6.33) 


RHS 
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KAr_Ar 

rAr) 

vx/' 

^7+1 


iAz  l 


4/cq/iq  \Ar+Ar  Ar_Ar 


iAz  /  2  1 


-  k%n  +  k20nl 


4koriQ  \Ar+Ar  rAr 


(6.34) 


This  is  commonly  referred  to  as  a  tridiagonal  matrix  system  of  equations  because  it  can  be  represented  in  matrix  form 
shown  in  section  6.5.1. 


6.3.3  Boundary  Conditions 

Equation  6.15  is  undefined  at  r  —  0.  To  circumvent  this,  axial  equations  can  be  added  to  the  system  of  equations  by 
applying  appropriate  boundary  conditions  to  the  problem.  Applying  the  boundary  conditions  from  equation  6.22,  S a 
(the  A  subscript  denotes  the  axial  solution)  becomes 


Sa  = 


2k 


2d 2 

dr2 


■  kl  \n2(r,z)  -  nl\ 


(6.35) 


Using  equation  6.30  and  the  new  operator  SA,  the  numerical  equation  at  the  axial  boundary  becomes  the  following: 


LHS 


t; 


;+i 


;+i 


iAz 

4^o«o  \Ar+ 
iAz  I  4 
4k0n0  \Ar+ 
iAz  I  2 


■  k^n  +  k^riQ 


4&o«o  \Ar+ 


RHS 


(6.36) 


up/ 


T' 


T", 


iAz  /  2 


4kotio  \A r2 
iAz  I  4 


iAz  l  2 


(6.37) 


4kofio  \Arr 


Applying  the  relationship  for  the  “phantom  point”,  4"_p  from  equation  6.21,  Equations  6.36  and  6.37  become  as  shown 
below: 


RHS 


T' 


% 


iAz 
4k0n0 
iAz  (  4 


7T  -  k0n  ■ 
Ar; 


•  klnl 


4koiio  \Ar+ 


(6.38) 
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LHS 


;+i 


t; 


'?'+1 


1  + 


iAz 


4-koiio  \A r\ , 

The  system  of  equations  created  can  be  represented  in  the  matrix  form  shown  in  section  6.5.1. 


(6.39) 


6.4  Pade  Method 


The  Pade  method  starts  with  equation  6.15  as  well.  In  this  case,  a  paraxial  solution  is  not  assumed  and  the  second 
derivative  in  the  axial  direction  is  left  intact.  Equation  6.15  is  rearranged  into  the  recursive  form  shown  in  equation 
6.40  to  become 


dz 


iP 

2k0n0 


1  + 


-T* 


2kon0  dz 


(6.40) 


where 


P  = 


d2  1  d  u,  q 

d?  +  ~rd~r+k 


(6.41) 


6.4.1  Pade  Approximate 

The  recursive  nature  of  this  equation  along  with  the  help  of  [1]  yields  the  Pade(l,l)  form  solution  in  Equation  6.42: 


d'V 

dz 


iP 

2k0n0  yjy 


1  + 


(6.42) 


Higher-order  Pade  approximates  are  available  such  as  the  Pade(3,3).  The  Pade(l,l)  derivation  will  result  in  solving 
a  tridiagonal  matrix  system  of  equations.  The  higher-order  approximates  result  in  solving  more  complex  systems  of 
equations  such  as  a  pentadiagonal  matrix  system  of  equations. 


Using  a  finite  difference  approximation  to  the  derivative  with  respect  to  z  on  the  left  hand  side  and  taking  an  aver¬ 
age  of  the  *¥  on  the  right  hand  side: 


VJ//+1  _  \J// 

A  z 


2^  yf+1  +  yf 

l  +  J-2  2 

4k~0n0 


By  rearranging  equation  6.43,  equation  6.44  can  be  obtained: 


(6.43) 


PT"'+1  +  iknoAzPY*1  +  4klnl'¥i+1  =  PT  -  ikn0AzPx¥l  +  4^1^  (6.44) 

Substituting  the  r  derivatives  and  constants  contained  in  the  P  operator.  Equations  6.45  and  6.46  can  be  found,  equa¬ 
tion  6.45  is  the  left  side  of  equation  6.44  and  equation  6.46  is  the  right  side  of  equation  6.44: 


LHS 
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dlxYi+l  1 d'Fi+1 

dr 2  r  dr 

d2yVM 


r/i+l 


+  ikn()Az 


dr2 


kl  [»2  -  nl\  'V' 

ki[n2-nl\' V- 


l  W" 

H - — h 

r  dr 


(6.45) 


+  4fcg«o'I'!+1 


RHS 


52vP'  1W  ,  2  r  2  21  uji 

+  -^+k«\n  -'i«r 


dr2 

-  ikiiyAz 


r  dr 

d2'V{  1  W  .  2  r  2  2  I  u/, 

77)7  +  r  _  ”°l  ^ 


<9r2 


(6.46) 


+  4/t2/i2vP'+l 


Applying  finite  difference  approximations  from  Equations  6.18  and  6.19  for  the  transverse  derivatives.  Equations  6.45 
and  6.46  become  as  shown: 


LHS 
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1 
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r 
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VI/i+ 1 

7+1 


VI/l+1 

■  /-I 


A  r 


Ar+ 
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1 

+  - 
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VI/i+1 
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A  r 


kl\n2-nl]%+1 
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vi/1+1 
'  7-1 


(6.47) 


+  4/t5«5vf"+1 
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(6.48) 


+  4/t,7i2vf"' 


By  rearranging  and  multiplying  by  A r.  Equations  6.47  and  6.48  become  the  equations  that  follow: 


LHS 


9YT//+1 

7- 1 

VI/i+1 

7-1 

2^o«oAzTi;_l|  ;l()n0Az'f"+11 

Ar_ 

r 

Ar_  r 

2vF71 

2'F'+1 

■  +  ^A/-[/r  «o|r;' 

Ar+ 

Ar_ 

2/*o«oA7F+' 
A  r+ 

+  4/fc2«2Ar'f''+l 


2ik0n0Az'¥iJ+l 

aTI 


ik0n0AzklAr  -  775]  'P^+1 
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- —  +  +  - —  +  - — 

A  r+  r  A  r+  r 


(6.49) 
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R I  IS 


24" 

7-1 

'Pf-i 

2ik0nQAzxi“i_ ,  ikonoAz'V1-  , 

Ar_ 

r 

Ar_  '  r 

2vf'(fl 

2xT'+1 

+  fc2Ar[«2-rc2]'P'+1 

Ar+ 

Ar_ 

UkonoAz'V)  UkonoAz^)  ,  r  ,  0l  . 

+ - — - -  + - — - -  -  ikonoAzk^Ar \>r  -  Tl 

+  4kln20Arx¥ij 

2^‘j+l  '¥ij+l  2ik0n0AzTj+[  ik(mAz^)+] 

A  r+  r  Ar+  r 

The  solutions  at  every  r  position  are  then  factored  out  of  Equations  6.47  and  6.50: 

LHS 


RHS 


2  1 
Ar_  r 

2  2 
Ar+  A  r- 


r/i+l 


T - |(1  +  ik0n0Az)%U 

k^Ar  yr  -  «q]  1(1+  ikon^Az)  +  Ak^n^Ar 


r; 


+  |  +  ^  |  (1  +  ik0n0Az)  T"+1 


j+ 1 


Ar_ 


(i  -/MoAz)r. ! 


- +  k^Ar  [n2  -  «o]j  (1  -  ikan0Az)  +  Ak^n^Ar 

[£r  +  ^(1-iMoAz)^.! 


6.4.2  System  of  Equations 

The  system  of  equations  in  6.51  and  6.52  can  be  rewritten  as 


a(-)/3(+y V‘;\  +  [ y/3(+)  +  A k20n20Ar\  T'f  +  a(+)JS(+)'P'++11 
=  «(-)/!(-)'?'_,  +  [7p(-)  +  Ak20n20Ar]  T'  +  rr(+)y8(-)'T'+1 

where 


(6.50) 


(6.51) 


(6.52) 


(6.53) 


a(±) 


2  1 

-  ±  — 

A  r(±)  r 


P(±)  =  1  ±  ik0n0Az 


7=  ~ 


2 

aZ 
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6.4.3  Boundary  Conditions 

At  r  —  0,  it  is  again  necessary  to  use  the  boundary  conditions.  Applying  equation  6.22,  the  operator  PA  becomes 


Pa  = 


2(f  if  2,  ,  9 1 

—  +lq\n-(r,z)-n-0  \ 


Applying  equation  6.4.3  makes  equation  6.44  become: 


PA T'+1  +  ikn0AzPA'Vi+'  +  4^^'Pi+1  =  Pa'V'  -  ikn0AzPA'i’i  +  4*^Ti  .  (6.54) 

Following  the  same  steps  as  for  the  rest  of  the  grid  points,  the  numerical  equation  at  the  r  —  0  boundary  becomes 


where 


[y/3(+)  +  4 k20n20Ar\  T'+1  +  a(+)P(+)%+1 
=  +  [y/5(-)  +  4k~n20Ar\  %  +  cx(+)P(-fV\ 


a(±) 


8 

Ar(±) 


fi(±)  -  1  ±  ik0n0Az 


7  ~ - t-  kz,Ar  \n 2  -  rix  1  . 

r  A  r+  0  L  °J 

Again  this  system  of  equations  can  be  written  in  the  matrix  form  found  in  section  6.5.1. 


6.5  Solving  the  Tridiagonal  Matrix  Equation 

6.5.1  Tridiagonal  Matrix  Equation 

The  system  of  equations  found  in  the  Caley  method  and  the  system  of  equations  found  in  the  Pade  method  can  be 
represented  in  the  following  tridiagonal  matrix  form.  The  matrix  values  are  the  terms  multiplied  by  the  'F  at  each  grid 
point.  These  matrix  values  are  different  for  each  method: 
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M- 1 


This  matrix  can  be  represented  in  the  matrix-vector  notation  as 


AT'  =  IBT'+1  . 


(6.56) 
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6.5.2  Explicit  Solution 

The  left  side  of  equation  6.56  is  made  of  known  information.  That  side  of  the  equation  can  be  solved  explicitly.  The 
computer  directly  multiplies  all  the  matrix  elements  with  the  known  solution  and  equation  6.56  becomes 

O  =  B'F,+1  .  (6.57) 

Where  ®  is  the  array  result  of  explicitly  solving  the  left  side  of  equation  6.56. 

6.5.3  Thomas  Algorithm 

Equation  6.57  can  not  be  solved  explicitly.  The  inverse  of  the  matrix  B  could  be  taken  and  then  explicitly  multiplied 
with  array  components  in  ®  to  find  VF,  but  this  is  computationally  intensive.  Instead,  the  Thomas  algorithm  [9]  is 
utilized.  The  Thomas  algorithm  is  a  fast,  efficient  way  of  solving  this  equation  without  finding  an  inverse.  ’’Numerical 
Recipes”  contains  an  algorithm  called  tridiag  that  was  used  for  this  step. 

6.6  Propagation  using  Hankel  Transforms 

6.6.1  Bessel  Functions 

Over  an  infinite  range,  [0,  °o],  the  Bessel  functions  form  a  complete  set.  Any  arbitrary  function  can  therefore  be 
represented  as  a  linear  combination  of  Bessel  functions.  The  set  of  Bessel  functions  can  be  written  as  abstract  ket 
vectors,  \JP.V),  of  the  vector  space  whose  inner  product  is  defined  as1: 

W)  =  2«f  f  (r) g  (r)  rdr  (6.58) 

Jo 

The  set  |  Jpy)  span  the  inner  product  space,  or  form  a  complete  basis.  The  index  p  denotes  the  Bessel  function  order, 
and  v  is  a  continuous  index.  Casting  this  into  the  coordinate  basis  gives 

(r\JP,v)  ->  JP  (27 zvr)  . 

Over  a  finite  range  [0,  /?],  the  continuous  index  is  quantized,  (v  — >  n).  The  complete  set  over  the  finite  range  becomes 

(r\JP,n)  ->  JP(^Tr)  n  =  1,2,...  , 

where  aPt„  denotes  the  nth  index  of  the  pth  order  Bessel  function.  The  Bessel  functions  are  convenient  to  work  with 
because  they  are  the  eigenfunctions  of  the  radial  component  of  the  cylindrical  Laplace  operator: 


or,  written  abstractly. . . 

ct' 

ro(r)|/p,n>  =  —j^\Jp,n) 

A  few  properties  of  the  Bessel  functions  will  be  useful  in  the  following  derivation.  Note  that  in  the  notation  used  here, 
OgpU P+\,v=i)  corresponds  to  Jp+\(ap^ 

Closure  (Continuous): 


(JpMpy)  =  -S(y-V) 

1  The  271  here  is  included  with  the  cylindrical  coordinates  in  mind. 


(6.59) 
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Orthogonality  (Finite): 


Normalization  (Finite): 


(Jp,n\Jp,ni)  —  0 


(J p,n\J p,n ) 


(6.60) 


(6.61) 


6.6.2  Projections 

The  projection  operators  of  the  Bessel  basis  and  the  coordinate  basis  are: 


P(J”) 

P(r) 


p,y)<J p,v\ 
p,v\J p,v) 

k)<r| 

(r\r) 


A  vector  in  the  space  can  be  represented  by  it’s  components  in  the  Bessel,  or  coordinate  basis  by: 


Is) 

Is) 


Note  here  that  over  the  infinite  range,  we  have 


P(4)|  g)dv 


P(r'|  g)dr 


<r|r> 


{J p,v\J p,v ') 


S(r-r') 

2nr 

6(v  —  v') 
2nv 


from  the  Closure  equation  (equation  6.59)  and  the  definition  of  the  inner  product  (equation  6.58),  and  the  projections 
of  g  can  be  written 


/-*oo 

Is)  =  2?r  I  \Jp,v)(Jp,v\g)vdv 
Jo 

Is)  =  27T  f 

Jo 


\r){r\g)rdr  . 


Over  a  finite  range,  the  continuous  index,  v,  on  the  Bessel  functions  becomes  discrete,  and  the  integral  will  go  to  a 
sum.  The  coordinate  basis  vectors  are  still  indexed  by  a  continuous  variable. 


Is)  =  2p(7")ls) 

n=  1 

rR 

Is)  =  IP(,)ls> 

Jo 

Again,  from  the  inner  product  definition  and  normalization  integral  (equation  6.61),  these  can  be  written 


Is) 


I J p-¥  1  ,V—  1 ) 


2 


Jn  IS) 


Is)  =  2 7T  f  \r){r\g)rdr  . 
Jo 
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From  these  representations,  it  is  obvious  that  the  sum/integral  over  projection  operators  is  the  identity  matrix. 

Two  approximations  will  be  used  concerning  the  sum  of  the  projection  operators.  The  hist  is  that  of  the  coordinate 
basis  projection  operators  over  a  finite  range.  The  operators  must  be  integrated  over  to  give  the  identity  matrix,  but 
can  be  approximated  by  a  sum: 


k)<r| 

(r\r) 


m= 0 


1^* m)(T ml 
ml F m) 


(6.62) 


Next,  we  note  that  if  the  index,  v,  has  a  finite  bound,  V,  then  the  infinite  limit  of  the  integral  can  be  replaced  with  V: 


\J p,v)\J p,v\ 
(j p,v  \J p,v ) 


si* 


\J p,v)^ p,v\ 
p,v  \J p,v ) 


These  will  be  used  as  approximations  to  the  identity  matrix. 


(6.63) 


6.6.3  The  Quasi-Discrete  Hankel  Transform 

The  beam  profile  is  assumed  to  be  an  azimuthally  symmetric,  separable,  function: 

'V  (r,  e,  z)  -»  'V  ( r ,  z)  -»  T'(,')  (r)  Tb)  (z) 

The  electric  held  distribution,  and  again  the  Bessel  functions,  can  be  written  as  abstract  ket  vectors:  2 


W 

J p,n  ) 

Casting  these  vectors  onto  the  coordinate  basis  retrieves  the  functional  representation 


(r,  z|'T)  -»  T(r,z) 

(r\Jp,n)  ^  Jp^r}  . 

If  R  is  chosen  to  be  larger  than  the  beam  radius,  the  radial  electric  held  distribution  can  be  written  as 


(r|'T)  =  <r|  2 

n=  1 


1*^ p,n}{J p,n\ 
(J p,n\J p,n ) 


1^), 


where  the  sum  of  the  projection  operators  has  been  inserted.  Slight  rearranging  of  terms  gives 


<r|'T>  = 

n=  1 


p.ii) 

p,n\J p,n) 


<jP*m  ■ 


(6.64) 


(6.65) 


It  is  apparent  that  the  sum  is  a  Bessel  series  expansion  of  the  electric  held  distribution.  The  inner  product  (,//v,|vI') 
is  the  electric  held  distribution  represented  in  the  Bessel  basis,  and  is  the  Hankel  Transform  of  the  original  electric 
held  distribution  represented  in  the  coordinate  basis,  (r|vF).  Similarly,  the  Hankel  Transform  can  be  expanded  in  the 
coordinate  basis: 


rR  (Jp,n\r)(r PF> 

Jo  <r|r> 


Using  equation  6.62  to  replace  the  integration  with  a  sum  and  rearranging  terms  gives 


m=  1 


p,n\rm) 
m\f m) 


<rmm  ■ 


“The  Bessel  functions  are  assumed  to  be  functions  of  the  r  coordinate 


(6.66) 
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Over  a  finite  range,  equation  6.66  is  an  approximation  of  the  Hankel  Transform,  while  equation  6.65  is  the  Inverse 
Hankel  Transform.  A  means  of  numerically  evaluating  these  transforms  can  be  reached  through  a  symmetric  deriva¬ 
tion. 


With  a  little  algebra,  the  two  representation  can  be  written  in  a  symmetric  form. 


RirJV) 


z 


2(7" m\J n  n) 


v(jP,nm 


|<^I4+1,V=1)|  tt  2nRV\{a^\Jp+hv=,)\\{a^\Jp+hv=l)\  |<Vi,v=il^)| 


|<^IVl,v=l>| 


p,n\f m) 


R{rmm 


5  2nVR  p+\,v=\)\  |('a?l,Vt-i,v=i)|  |<Vi,v=tl^)| 
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We  define  two  vectors  and  a  matrix. 


(6.67) 

(6.68) 

(6.69) 

and  the  two  transformations  can  be  written  as  matrix  multiplication: 

<DW  =  Tm,n<D(/) 

4>(y)  =  Tn,m<D(r) 


|(^IVl.v=l>| 

(y)  _  njp,  m 
"  |<lfl-W=i>| 

m\J p,n ) 

m’"  ”  2^y|<^|7p+i,v=1)||<^|7p+i,v=1)| 


6.6.4  The  Propagator 

Once  the  radial  electric  field  distribution  is  known  in  terms  of  it’s  Bessel  Function  components,  a  method  referred  to 
as  the  Split-Step  Method  is  used  to  solve  the  Wave  Equation.  Equation  6.15  is  split  into  two  equations. 


d2xV  .  d'V  d2^  1  d'V 

~2lk°n°~dl~dr2  + 


+  kl[n2(r,z ) 


d 2xf'  ff¥  32vT  1  d'V 

—  -2ik0n0—  +  ^  +  -  — 

d2xV  W  ,r  , 

~^Y  -  2ik0n0—  +  k~Q  [ n  (r,  z)  -  «0J  ^ 


=  0 


i 


=  0 
=  0 


(6.70) 

(6.71) 

(6.72) 


and  the  two  are  solved  separately.  Notice  that  the  equation  6.70  is  just  the  Wave  Equation  in  a  homogeneous  media, 
and  equation  6.71  contains  the  inhomogeneous  refractive  index  terms.  The  strategy  will  be  to  solve  equation  6.70,  and 
step  forward  in  the  z  direction  by  some  Az.  Equation  6.71  can  be  interpreted  as  a  phase  correction,  and  will  be  used  to 
adjust  the  phase  of  the  wavefront,  based  on  the  refractive  index  profile  it  passed  through  during  the  Az  step.  Expanding 
the  radial  electric  field  distribution  in  terms  of  Bessel  Functions,  equation  6.70  can  be  rewritten. 
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With  a  little  manipulation,  the  r  differential  operator  on  the  right-hand  side  can  be  replaced  by  eigenvalues. 

d2 


q2  q 

— 7  -  2ik0n0— 
dz~  dz 


V  (z)  V  akJp  +  T*  (z)  Y  ak 

ti  K  R  ’  ti 


1  d_ 

dr2  r  dr 


=  0 


I 

Jfc=l 


C-fr)  = 


The  Bessel  Functions  are  linearly  independent,  so  each  index  i  must  separately  equal  zero, 

\2 


92  ry,  8 


(6.73) 
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Now  Equations  6.70  and  6.73  have  the  same  form.  In  general,  the  solution  to  these  differential  equations  is: 


C\2  O 

-Q2  -  2ik0n0  -7-  f  (z)g(r)  —  K2f  (z)  g  (r)  =  0 

/  (z)  g{r)=f  (z0)  g  (r)  V(W^ 

If  the  function  is  known  at  some  z  coordinate,  the  propagator  gives  the  function  at  all  other  z  coordinates.  Applying 
this  to  Equations  6.70  and  6.73,  an  incident  beam  profile  can  be  propagated  a  distance  A z  into  the  medium. 


r  (Az)T'  (r)  = 

4*  (Az)  4»  (r)  =  'Y  (0)  4''  (r)  e^hn(r,z) 


— r)  e~ik°"0eiAz  V  (<r°"o)2'(  ¥ )’ 


(6.74) 

(6.75) 


First  the  wavefront  is  propagated  a  distance  A z  with  equation  6.74,  then  Equation  6.75  is  used  to  adjust  the  wavefront 
phase  to  account  for  the  refractive  index  profile  in  the  slice  Az. 

In  terms  of  the  vectors  defined  in  6.67  and  6.68,  this  can  be  written: 


(7  +  a  z)  =  ©(-0  (z)  e~ikon°  eiAz  V (i°"o)2_(  t2? 

©£?  (z  +  Az)  =  (z)  e-'10"0  eiAzkon(r'z) 
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Chapter  7 

Monte  Carlo  Scattering 


Source  terms  within  highly  scattering  media  are  computed  using  a  Monte  Carlo  simulation.  The  Monte  Carlo  approach 
to  photon  transport  uses  a  statistical  method  to  compute  the  photon  distribution  within  a  tissue  that  may  be  scattering 
the  photons.  The  implementation  in  the  BTEC  model  utilizes  work  by  Wang  and  Jacques  [12]. 


7.1  Theory 

As  mentioned,  the  Monte  Carlo  method  for  modeling  photon  transport  is  based  on  statistical  distributions  and  random 
numbers.  Photons  are  essentially  modeled  as  ballistic  objects.  A  packet,  or  bundle,  of  photons  with  weight,  w,  is 
launched  into  a  tissue  and  propagated  like  ping-pong  balls,  randomly  interacting  with  the  tissue. 

Photon  propagation  is  done  on  a  3-dimensional,  Cartesian  coordinate  system.  A  photon  packet  has  a  position  vector, 
given  by  the  coordinate  vector,  r,  in  the  global  coordinate  system,  and  a  unit  direction  vector,  n,  that  specifies  the 
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direction  the  packet  is  currently  traveling.  The  packet  propagates  a  distance,  s,  along  the  direction  vector  from  r,  to  rf, 
where  rf  is  given  by 

rf  =  ri  +  sn  .  (7.1) 

After  each  propagation,  the  packet  interacts  with  the  tissue,  where  some  of  the  photons  are  absorbed,  corresponding 
to  a  decrease  in  the  packet  weight,  and  the  rest  are  scattered.  The  packet  is  propagated  in  the  new  direction  and  the 
process  is  repeated  until  all  photons  in  the  packet  are  absorbed. 

The  photon  absorption  is  recorded  on  a  2-dimensional  cylindrical  grid.  At  each  interaction  site,  the  closest  grid  point 
is  determined  from  the  r  =  fx 2  +  y 2  and  z  coordinates  of  the  photon  packet,  and  the  absorbed  packet  weight  is  tallied 
at  the  grid  location.  If  the  packet  is  determined  to  be  outside  the  boundaries  of  the  grid,  then  it  is  killed,  and  a  new 
packet  is  launched.  Any  remaining  weight  is  added  to  one  of  three  tallies,  indicating  which  external  boundary  the 
packet  crossed  (reflected,  transmitted,  or  lateral). 

As  more  photon  packets  are  launched  into  the  tissue,  a  photon  weight  distribution  is  compiled  on  the  grid,  representing 
the  number  of  photons  absorbed  at  (or  near)  each  grid  point  in  the  computational  space.  The  distribution  is  initially 
very  discontinuous,  but  smooths  out  as  more  photon  packets  are  propagated.  Once  photon  propagation  is  complete, 
this  distribution  is  transformed  into  a  photon  density  distribution,  differing  from  the  photon  weight  distribution  because 
of  the  cylindrical  coordinates  of  the  grid.  The  photon  density  distribution  is  normalized  to  the  power  of  the  simulated 
emitter  (minus  the  percentage  of  photon  weight  that  was  not  absorbed  on  the  grid).  The  resulting  distribution  gives  the 
source  term  for  the  heat  equation. 

The  amount  of  computation  required  for  the  Monte  Carlo  simulation  varies  greatly  with  the  model  configuration.  As  a 
general  rule  though,  the  amount  of  time  a  simulation  requires  increases  with  higher  scattering  coefficients,  and  lower 
absorption  coefficients.  This  happens  to  be  the  two  instances  in  which  the  model  is  needed  most,  so  the  simulation 
time  is  inherently  long.  Therefore,  photon  propagation  is  only  computed  when  necessary.  The  Monte  Carlo  routine 
only  runs  when  the  emitter  is  actually  turned  on  in  the  thermal  simulation.  After  the  source  term  is  computed,  it  is 
saved,  and  reused  for  future  time  steps  in  the  heat  model.  This  method  is  valid  as  long  as  there  are  no  changes  in  the 
thermal  model  that  will  affect  the  propagation  simulation  (such  as  temperature  dependent  absorption  and  scattering 
coefficients). 


7.1.1  Photon  Step  Size 


The  distance  a  photon  travels  between  interaction  sites  is  determined  by  sampling  the  probability  distribution  of  the 
photon’s  free  path,  s.  After  each  tissue  interaction,  the  photons  free-flight  step  size  is  determined  using  the  interaction 
coefficient  [12], 


~ln(Q 

Hr 


(7.2) 


where  f  is  a  random  number  between  0  and  1  that  is  uniformly  distributed.  The  interaction  coefficient  gives  the 
probability  interaction  per  unit  length,  and  is  related  to  the  absorption  and  scattering  coefficients  of  the  tissue: 


1 

Hr  =  - 7 - 

Ha  +  Hs 


(7.3) 


Note  that  the  term  - ln( %)  can  be  interpreted  as  a  dimensionless  free-flight  distance.  The  actual  step  size  is  computed 
by  dividing  this  dimensionless  quantity  by  the  interaction  coefficient  of  the  tissue.  If  the  photon  propagates  across 
multiple  tissues,  the  dimensionless  step  size  must  be  scaled  by  the  interaction  coefficients  of  each  tissue  propagated 
through.  The  actual  free-flight  distances  through  each  tissue  are  related  to  the  dimensionless  step  size  by 


Y  HnSi  =  ~ln(Q  ,  (7.4) 


where  /j,i  is  the  interaction  coefficient  of  the  /th  tissue,  and  .v,-  is  the  distance  traveled  in  that  tissue.  Boundary  interac¬ 
tions  will  be  discussed  in  section  7.1.3 


50 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


7.1.2  Tissue  Interaction:  Absorption  and  Scattering 

After  a  step  size  is  chosen,  the  photon  packet  is  propagated  that  distance  along  the  direction  vector.  At  the  new 
position,  part  of  the  packet’s  weight  is  absorbed.  The  absorbed  weight  is  computed  as 

dw  -  —w  .  (7.5) 

Ft 

The  weight  of  the  photon  packet  is  decreased  by  this  amount,  while  the  total  weight  at  the  nearest  grid  point  is 
increased. 


w  -  w  -  dw  (7.6) 

A[i][j]  =  A[i][j]  +  dw  (7.7) 

Here  A  has  been  used  to  denote  the  array  corresponding  to  the  cylindrical  grid  storing  the  photon  weight  distribution. 
Next  the  packet  is  scattered,  which  just  corresponds  to  a  change  in  the  packet’s  direction  unit  vector.  Scattering 
involves  two  angles:  a  deflection  angle,  ©  — »  [0,  n],  and  a  rotation  angle,  <D  — >  [0, 2n],  The  anisotropy  of  the  tissue,  g. 


Figure  7.2:  Photon’s  direction  n  is 
(Photon  not  to  scale) 

is  a  measure  of  the  expected  value 


deflected  by  an  angle  ©  from  the  propagation  axis,  and  rotated  about  the  axis  by  (£>. 


of  cos©.  The  probability  distribution  of  cos 0  is  given  by 

i-r 


P ( cos  ©)  = 


2(1  +  g2  -  2  g  cos  ©)3/2 


(7.8) 
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so  for  each  scattering  event,  cos ©  can  be  randomly  chosen  using  a  uniformly  distributed  random  number,  between 
0  and  1 : 


cos  0 


2« 


i-r 

l-g+2gf 


if  8*  0 


(7.9) 


2f-  1 


if  8  =  ® 


The  second  angle,  <D>,  has  a  uniform  probability  of  being  between  0  and  2 n.  The  photon  is  just  as  likely  to  scatter  left, 
as  it  is  to  scatter  right: 

i//  =  2ni;  (7.10) 


A  new  step  size  is  chosen,  and  the  packet  is  propagated  along  the  new  direction. 


7.1.3  Boundary  Interaction:  Reflection  and  Refraction 

The  photon  packet  has  a  possibility  of  intersecting  a  tissue  boundary  during  a  propagation  step.  If  so,  the  packet 
should  not  be  propagated  across  the  boundary  transparently  for  two  reasons.  One  being  that  the  step  size  chosen  is 
only  correct  for  the  tissue  at  the  packet’s  initial  position.  If  the  packet  propagates  into  another  tissue,  the  step  size 
must  be  rescaled  for  the  new  tissue.  Secondly,  if  the  refractive  index  of  the  two  tissues  differs,  the  packet  will  either 
be  reflected  off  the  boundary,  or  refracted  across  the  boundary. 

The  probability  that  the  packet  is  reflected  at  a  given  incident  angle  is  given  by  the  reflection  coefficient: 


R(0i)  = 


1 

2 


sirr(9j  -  6 ,) 
sin2(0j  +  0t) 


tarf(9i  -  9,) 
tan2(9i  +  9,) 


(7.11) 


When  a  packet  reaches  a  boundary,  a  random  number,  uniformly  distributed  between  0  and  1,  is  chosen.  The 
incident  angle,  9,,  is  then  computed.  If  £  <  R  (0,),  the  packet  is  reflected,  otherwise  the  packet  is  transmitted. 

With  tissue  layers  being  perpendicular  to  the  z-axis,  reflection  simply  corresponds  to  a  change  in  sign  of  the  z- 
component  of  the  photon’s  direction  vector,  and  refraction  is  given  by  Snell’s  Law: 


riisin(9\ )  —  msin  (tL) 


Care  must  be  taken  when  9  approaches  0.  If  the  photon  travels  some  distance,  d\,  before  intersecting  the  tissue 
boundary,  then  the  distance  left  to  travel  is  t  =  Si  —  d\.  If  the  photon  is  reflected,  propagation  continues  in  the  new 
direction,  with  the  new  step  size,  t.  If  the  photon  is  refracted  across  the  tissue  boundary,  the  travel  distance  must  be 
scaled  for  the  new  tissue.  Initially,  a  dimensionless  step  size  is  chosen,  and  scaled  by  the  interaction  coefficient  of  the 
tissue  to  give  an  actual  step  size  (equation  7.2).  The  dimensionless  quantity  for  the  distance  left  to  travel  is  given  by: 

s'  =  fi,i  (si  -  d i) 

This  quantity  is  scaled  by  the  interaction  coefficient  of  the  second  tissue  to  give  the  new  step  size: 

f  Bn  ,  ,  . 

s2  =  —  =  —  (si  -d\) 

Rt2  Bl2 

After  the  boundary  interaction,  propagation  continues.  If  other  boundaries  are  hit,  the  process  is  repeated. 


7.1.4  Packet  Termination 

Packets  are  terminated  when  their  weight  drops  below  a  threshold.  The  threshold  is  set  to  be  0.01%  of  the  initial  packet 
weight.  However,  the  packet  cannot  simply  be  terminated  once  the  weight  drops  below  threshold.  To  conserve  energy, 
a  method  called  Russian  Roulette  [12]  is  used  to  terminate  the  packet.  If  the  packet  weight  is  below  threshold  after  an 
absorption  event,  it  enters  a  game  of  roulette.  The  packet  has  a  1  in  10  chance  of  surviving  the  game.  If  it  survives,  it 
propagates  another  step,  but  will  have  to  play  roulette  again.  If  it  does  not  survive  the  roulette,  it  is  terminated  and  a 
new  packet  is  launched. 


52 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


7.1.5  Calculating  the  Source  Term 

After  a  complete  Monte  Carlo  simulation,  we  are  left  with  a  grid  that  contains  the  photon  count  at  each  point.  This 
photon  count  needs  to  be  turned  into  a  power  density  corresponding  to  the  emitter  power.  The  grid  is  a  2-dimensional 
slice  from  a  3-dimensional  cylinder.  Therefore,  the  photon  counts  at  each  grid  point  represent  the  totals  for  a  ring  in 
the  cylinder.  At  a  specific  grid  point,  A  [/][/],  the  volume  of  the  corresponding  ring  is 


V  -  n 


Az+  +  Az 

r]  2  J 

2 

The  photon  density  on  the  2-dimensional  cylindrical  grid  is  computed  from  the  photon  totals  as 


mij] 


A[i][j] 
V[i][j ]  ’ 


where  V[t][j]  is  the  corresponding  volume  element  for  the  point  A  [z][y].  A[i][y]  contains  the  photon  density,  which 
is  directly  proportional  to  the  power  density.  During  photon  propagation,  four  quantities  are  tracked:  total  weight 
launched,  total  weight  reflected,  total  weight  transmitted  and  total  weight  lost  laterally.  The  ratio. 


wt  -  (wr  +  w,  +  W/) 

a  =  - 

w, 


is  the  ratio  of  power  deposited  in  the  tissue  to  power  input.  To  compute  the  power  density,  the  photon  density  is 
normalized  to  the  amount  of  power  deposited  in  the  tissue  by  the  emitter. 


aPe  —  2  n 


A  (r,z)  r2drdz  ■ 


Integration  is  taken  over  the  limits  of  the  computational  space.  After  normalization,  the  grid,  A[;][j],  contains  the 
necessary  source  term  for  the  heat  equation. 


0  0.05  0.1  0.15  0.2 

z  (cm) 
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z  (cm) 


Figure  7.3:  Two  computed  source  terms,  with  (left)  and  without  (right)  scattering. 
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Figure  7.4:  Scattering  ^progression  and  logic 
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Chapter  8 


Z-scan  Simulation 


8.1  Background  and  Geometry 

A  z-scan  experiment  is  commonly  used  to  observe  thermal  lensing.  In  a  z-scan  experiment,  a  laser  is  focused  and  the 
output  is  observed  through  a  circular  aperture.  A  cuvette  sample  is  placed  at  different  distances  on  each  side  of  the 
focal  point.  Water  will  be  used  for  this  discussion,  although  other  media  may  be  used.  The  requirement  is  that  the 
material’s  refractive  index  changes  with  temperature.  As  temperature  rise  in  the  sample  causes  the  beam  to  refocus 
because  of  the  refractive  index  change,  the  power  transmission  through  the  aperture  is  recorded.  This  is  done  with  the 
sample  at  many  locations  between  the  lens  and  the  aperture  along  the  beam  axis,  hence  ”z-scan.” 


Screen 


Figure  8.1  is  a  simplified  example  of  the  z-scan  experiment.  In  this  figure,  the  laser  beam  comes  in  from  the  left.  It 
passes  through  geometric  optics  and  then  through  the  water  sample.  The  beam  passes  through  a  fixed  circular  aperture 
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and  then  is  either  observed  on  a  screen  or  the  power  is  measured  with  a  detector.  Again,  the  water  sample  is  moved  to 
many  positions  between  the  lens  and  the  circular  aperture. 

The  BTEC  model  is  equipped  to  simulate  the  z-scan  experiment.  It  uses  the  Gaussian  Beam  Propagation  model 
discussed  in  chapter  5  to  calculate  the  electric  field  through  the  linear  optics  and  atmosphere  up  until  the  dispersive 
sample.  The  finite  difference  beam  propagation  (chapter  6)  model  and  the  two-dimensional  heat  model  (chapter  ??) 
are  used  together  to  propagate  the  laser  through  the  sample  and  compute  the  temperature  rise.  To  the  right  of  the 
sample  in  figure  8.1,  the  Kirchhoff  diffraction  integral  is  used  to  calculate  the  electric  field  intensity  at  all  points.  The 
model  can  show  how  the  electric  field  would  look  projected  on  a  screen  (see  the  top  diagram  of  figure  8.1)  and  can 
compute  percent  transmittance.  The  diffraction  code  numerically  approximates  the  percent  transmittance  based  on  the 
electric  field  distribution  at  the  aperture  plane  and  the  aperture  diameter  defined  in  the  configuration  file. 


8.2  Derivation  of  Diffraction  Integral 

The  electric  field  diffraction  pattern  on  a  plane  a  distance  d  from  from  the  imaging  plane  is 

E(R,  0)  =  —e~ik°RD(0) 

Ad 

R  =  -yjq2  +  d2 

where  q  is  the  radial  distance  on  the  imaging  plane. 

Integrate  the  function  over  all  radial  input  points  and  then  integrate  again  over  all  q  on  the  imaging  surface: 

E  =  E(r) 


(8.1) 

(8.2) 


0(0)  = 


(8.3) 


where  R  is  not  a  function  of  r. 

Express  in  terms  of  E,  noting  that  it  is  complex: 

in  rq  a  rnP~l  l2nr\ 

E(R,0)lmag  =  Yde~'k°R  J(}  Rq  J0  Irn(E)J0\  —  \rdr  (8.4) 

in  it  r  Cq  0  rnP~l  (271  r\ 

E(R,0)rea/= —e-,knR  J  pdq  J  Re(E)J0l—rdr,  (8.5) 

where  Jq()  is  a  first-order  Bessel  function  found  by  a  Numerical  Recipes  [9]  algorithm  called  BessJo. 

In  order  to  get  power,  E  is  multiplied  by  its  complex  conjugate  and  integrated  over  all  radial  points.  The  variable, 
E,  here  is  still  a  function  of  the  radius  and  theta. 

E  x  E*  =  | E(R,  0freal\  +  | E(R,  CJ  (8-6) 

So  from  the  expression: 


I(R,  0)  =  A  f  EAR,  6)E(R,  0)rodr  , 

Jo 

A  is  a  normalization  constant  with  the  extra  r  for  the  differential  in  cylindrical  coordinates: 


(8.7) 
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where  Pj  is  the  total  power. 


W,9)m 


E*Er0dr  =  PT 


Integrate  in  parallel  over  all  points  up  to  the  limiting  aperture: 


m,9)ap 

where  Pa  is  the  power  with  the  limiting  aperture. 


E*Er0dr  =  PA 


Divide  to  get  fractional  intensity: 


Pa 

Pt 


=  D(rap ) 


This  is  the  fractional  intensity  (W  or  J)  through  the  measurement  aperture. 


(8.8) 


(8.9) 


(8.10) 


8.3  Results 


A  comparison  between  a  z-scan  experiment  and  the  BTEC  model  is  shown  in  figure  8.2. 

Z-Scan  2.0  mm  Sample  (Water)  100ms 

1.4  r 


0  20  40  60 

Z-Scan  Position  (cm) 


80 


Figure  8.2:  Comparison  between  model  and  experiment 


The  experiment  was  preformed  with  a  1315-nm  laser.  The  water  sample  was  0.2  cm  thick.  This  is  initial  comparisons 
with  some  parameter  adjustment  to  estimated  material  parameters.  Each  data  point  was  taken  with  the  sample  at  a 
different  position  along  the  optical  axis.  At  that  position  the  transmittance  was  recorded  at  100  ms. 
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Chapter  9 


Damage 


9.1  Damage  Integral 

The  current  BTEC  heat-transfer  and  thermal  injury  model  produces  a  solution  to  the  2-D  heat  equation  that  provides 
an  approximate  time-dependent  solution  of  the  form,  T{z,  r,  t),  giving  the  temperature  distribution  in  cylindrical  coor¬ 
dinates  at  a  specific  time.  At  each  point  within  the  computational  space,  the  model  continually  updates  a  running  sum 
of  the  Arrhenius  integral  given  below: 


nfcr,-r)  =  AjT  (9.1, 

In  the  computational  space  z  and  r  are  indexed  by  (7,  j)  in  a  fixed  coordinate  space  such  that  the  values  may  be  indexed 
as 

<92> 

Currently,  this  equation  is  evaluated  employing  a  simple  trapezoid  method  which  evaluates  the  integrand  at  the  current 
and  previous  time  step,  maintaining  a  running  sum  of  the  integral  with  time.  For  a  given  coordinate,  the  i,  j  labeling  is 
ignored,  yielding 


O  (t  +  At)  —  Q  (r)  +  - 


exp  f 


-Ea 


RT(t  +  At) 


■  exp 


-Eg 

RT(t) 


(9.3) 


This  method  suffers  from  errors  as  the  integrand  varies  rapidly,  and  the  two  points  of  evaluation  are  weighted  equally, 
essentially  linearizing  an  exponential  function.  If  this  function  varies  such  that  the  equality  in  equation  9.4  is  not  a 
valid  approximation,  then  the  integral  is  not  accurate.  In  a  stretched  time-step  scheme  where  the  temperature  is  rising 
rapidly  at  the  end  of  a  pulse,  but  near  the  estimated  damage  threshold,  there  may  be  significant  error  as  that  time  step 
contains  nearly  all  of  the  running  summed  value  for  equation  9.3: 


exp 


RT(t) 


1  + 


(^)T 


(9.4) 


Obviously,  improved  methods  are  needed  to  accurately  depict  short-pulse  damage  thresholds,  where  temperature  rise 
is  significant  during  a  time  step. 


9.2  Threshold  Search 

In  the  current  BTEC  heat-transfer  and  thermal  injury  model,  an  optional  search  algorithm  can  be  used  to  determine  the 
threshold  of  an  effect.  The  search  algorithm  adjusts  the  power  of  an  emitter  until  the  threshold  criterion  is  met  to  within 
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some  user-specified  tolerance.  The  optional  search  can  determine  the  emitter  (or  combined  emitter)  power  required 
to  reach  a  user-specified  value  of  damage  threshold  or  a  user-specified  value  of  temperature  increase.  These  may  be 
global  criteria,  or  specified  for  a  certain  radial  extent.  For  example,  the  user  may  choose  to  evaluate  the  threshold 
power  required  to  obtain  a  damage  integral  value  of  1.0  anywhere  within  the  tissue,  or  at  a  specified  radius,  such  that 
a  lesion  size  may  be  used  as  the  threshold  determining  factor. 

The  current  search  algorithm  is  a  midpoint-search  method,  with  the  initial  search  range  (max,  min)  specified  as  a 
user-defined  input.  This  method  continually  bisects  the  search  range  until  the  threshold  value  is  reached.  In  practice, 
this  method  requires  about  6-12  iterations  in  which  the  BTEC  Thermal  Model  is  run  through  the  entire  simulation  time. 
The  bisection  method  was  selected  due  to  the  strongly  non-linear  nature  of  the  damage  integral,  for  which  Newton  or 
secant  methods  were  unstable. 

While  the  current  search  method  is  reliable,  the  currently  large  workload  for  the  model  running  many  high-resolution 
serial  jobs  on  the  AFRL/RHDO  clusters  is  significant.  This  is  creating  typical  1-2  day  time  for  a  job  to  reach  execution 
within  the  queues.  Any  reduction  in  search  time  would  be  significant  and  greatly  improve  turn-around  and  efficiency. 
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Chapter  10 

Stability  and  Accuracy 


The  question  of  stability  is  very  important  when  dealing  with  numerical  methods.  Much  of  the  discussion  of  stability 
only  applies  if  the  problem  space  is  one  homogeneous  material  with  sink  boundary  conditions.  The  cost  of  having  dif¬ 
ferent  layers  of  materials  or  more  complex  boundary  conditions  introduces  the  possibility  of  weakening  the  numerical 
strength  of  the  method.  The  most  important  thing  is  to  make  sure  that  the  methods  are  not  pushed  so  hard  that  they 
become  unstable.  Fortunately,  there  are  some  rules  that  should  keep  this  from  happening  in  most  cases. 


10.1  Different  Types  of  Errors 

Before  looking  at  the  method  it  is  important  to  know  about  the  two  main  types  of  errors.  The  total  error  will  be  the 
sum  of  these  two  errors. 

Round-off  Error 

The  first  type  of  error  is  round-off  error.  This  happens  because  of  a  computer’s  inability  to  exactly  store  floating  point 
numbers.  This  type  of  error  is  normally  very  small  and  does  not  cause  a  problem.  However,  it  has  been  seen  to  effect 
the  solution  if  the  solution  is  changing  very  rapidly.  That  is,  if  the  maximum  value  of  the  solution  changes  by  an  order 
of  magnitude  in  two  or  three  iterations  of  the  method.  If  the  solution  becomes  chaotic  and  evolves  very  rapidly,  it  is  a 
good  idea  to  check  if  it  was  rounding  error.  The  best  way  of  checking  this  is  to  change  the  data  type  and  see  if  where 
the  solution  becomes  chaotic  changes.  A  solution  that  started  to  look  chaotic  on  the  sixth  iteration  using  the  double 
data  type  should  take  more  iterations  to  see  the  chaotic  behavior  using  the  long  double  data  type.  Reducing  the  size  of 
the  time  step  should  stop  the  effect  of  round-off  error. 

Truncation  Error 

Truncation  error  comes  from  using  finite  difference  in  place  of  derivatives.  This  type  of  error  is  represented  using  the 
big  “O”  notation.  In  the  BTEC  code,  the  space  and  time  finite  differences  are  of  order  0(h2)  when  using  a  constant 
grid  spacing,  h.  The  one-dimensional  model  has  all  the  correct  finite  differences  in  space  to  handle  stretching  the  grid 
and  stay  0(h2).  Of  course,  the  h  gets  bigger  with  stretching  so  the  error  goes  up.  The  two-  dimensional  model  does 
not  use  a  second  derivative  finite  difference  that  handles  a  stretched  grid  so  its  truncation  error  will  not  stay  0(h 2)  as 
the  grid  is  stretched. 


10.2  Stability  of  the  Thomas  Algorithm 

Both  the  methods  used  in  the  BTEC  code  use  the  Thomas  algorithm  to  solve  a  tridiagonal  system  of  linear  equations. 
The  Thomas  algorithm  is  a  non-pivoting  elimination  method  that  is  very  efficient  for  solving  tridiagonal  systems.  A 
tridiagonal  system  is  of  the  form 
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Ax  =  b  (10.1) 

where  A  is  a  known  matrix  with  three  diagonals,  x  is  an  unknown  vector,  and  b  is  a  known  vector.  The  expanded 
tridiagonal  system  can  be  seen  in  equation  10.2. 


b\  ci  0 
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There  are  three  requirements  on  the  terms  in  matrix  A  to  keep  the  method  stable. 


(i)  a,  <  0,  bj  >  0,  and  c,  <  0 


(ii)  hi  >  -(a,+ 1  +  a- 1)  for  i  =  1,2  , . . . , A,  defining  co  =  a/v+i  =  0 

(iii)  b i  >  -( a,-  +  c,)  for  i  =  1,2 ,N,  defining  =  c#  =  0 


(10.2) 


For  a  proof  of  this,  please  see  [10].  It  will  be  shown  that  these  requirements  effect  boundary  conditions  and  put  a  limit 
on  the  how  much  properties  of  materials  can  change  from  layer  to  layer. 


10.3  Crank-Nicholson  Method 


The  Crank-Nicholson  Method  is  used  for  the  one-dimensional  model  in  the  BTEC  code.  The  Crank-Nicholson  Method 
has  an  accuracy  of  0(Ax2 ,  At2)  in  time  and  space.  It  is  also  said  to  be  unconditionally  stability  [7],  Unconditionally 
stable  means  there  is  no  limit  on  the  value  of  r  where 


At  K 

(A  x)2p  c 


(10.3) 


The  “no  limit”  means  that  r  can  be  as  large  as  desired  and  the  solution  will  not  blow  up.  This  does  not  mean  that  the 
solution  will  be  the  correct  solution  or  even  make  physical  sense.  To  insure  a  meaningful  solution,  it  has  been  shown 
that  there  is  a  bound  on  r  [10].  This  bound  is 


r  < 
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(10.4) 


where  L  is  the  length  of  the  model.  This  can  be  restated  as  a  limit  on  At  since  that  is  the  easiest  parameter  to  change. 


LAxp  c 
At  <  - — 

K7t 


(10.5) 


10.3.1  Layers 

The  previous  discussion  assumes  no  changes  in  material  properties  in  the  problem  space,  which  will  be  referred  to  as 
layers.  Nor  does  it  give  much  information  on  the  effect  of  boundary  conditions.  To  find  the  limits  of  the  method  with 
boundary  conditions  and  layers,  the  requirements  on  the  Thomas  algorithm  must  be  applied  to  the  method. 

By  applying  requirement  (i)  of  the  Thomas  algorithm  to  equation  2.25,  the  following  is  obtained  using  equations 
(2.13-2.15): 


2a-,-  >  6xKiAxi+ 


(10.6) 
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(10.7) 


2 Kt  >  -  6xKiAx/- 


2 piCi  SxKj[Axl+  -  A Xi-\  -  2k, 

~aT  >  A.v,  A.v,.  '  (10'8) 

The  only  way  any  of  these  equations  will  violate  the  requirements  of  the  Thomas  algorithm  is  if  6xk,  4-  0.  This  can  only 
happen  if  k,  is  changing,  which  occurs  on  the  boundary  of  where  two  layers  meet.  Equation  10.8  will  be  disregarded 
because  either  equation  10.6  or  equation  10.7  will  be  violated  first.  Equation  10.6  and  equation  10.7  provide  a  limit  on 
how  much  the  thermal  conductivity  can  change  going  from  one  layer  to  another.  They  can  be  simplified  by  assuming 
a  constant  grid  spacing  to 


4 Kj  >  KM  -  Kj_  | 


(10.9) 


and 


4/f,  >  Ki-\  -  Ki+l.  (10.10) 

This  means  that,  as  a  rule  of  thumb,  the  difference  in  thermal  conductivity  between  the  two  layers  needs  to  be  at  least 
less  than  four  times  the  smallest  of  the  two  for  the  Thomas  algorithm  and  thus  the  Crank-Nicholson  Method  to  be 
stable.  Running  the  Crank-Nicholson  Method  shows  that  it  becomes  unstable  when  this  rule  of  thumb  is  broken. 


10.3.2  Boundary  Conditions 

The  requirements  for  the  Thomas  algorithm  were  checked  for  boundary  conditions  and  it  was  found  that  the  boundary 
conditions  do  not  violate  the  requirements.  In  addition,  according  to  [10]  for  derivative  boundary  conditions  the  Crank- 
Nicholson  equations  are  unconditionally  stable.  For  the  convective  boundary  condition  a  large  number  of  solutions  for 
various  values  of  he  and  time  can  be  found  in  a  transient  temperature  chart  in  figure  2-13  of  [6].  These  solutions  can 
be  used  to  check  the  error. 


10.4  Peaceman-Rachford  Method 

The  Peaceman-Rachford  method  is  unconditionally  stable  and  has  an  accuracy  of  0(Az2 ,  Ar2 ,  At2)  [8].  The  same 
conditions  layer  interfaces  that  were  found  for  the  Crank-Nicholson  Method  can  be  found  for  the  Peaceman-Rachford. 
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Chapter  11 


Verification 

11.1  Introduction 


In  order  to  verify  the  BTEC  model,  several  validation  cases  were  considered.  Model  predictions  were  compared  to 
known  analytical  solutions  from  the  text.  In  some  cases,  predictions  were  also  compared  experimental  results,  data 
from  text,  or  with  other  implementations  of  the  specific  portion  of  the  model. 


11.2  Source  Term  Verification 


To  verify  the  calculation  of  the  source  term  used  in  the  2-D  heat  model,  the  BTEC  is  run  for  a  short  time  with  a  flat-top 
beam.  The  source  term  follows  Beer’s  Law.  The  beam  chosen  for  source  term  verification  is  0.1  cm  in  diameter  has 
100  W  of  power,  and  a  wavelength  of  1315-nm.  The  beam  was  set  to  a  zero  divergence.  The  beam  is  propagated 
through  water  for  100  /-is.  Absorption  of  a  1315-nm  beam  in  water  is  assumed  to  be  1.33  cm'1. 

The  initial  intensity  is  calculated  by 


(11.1) 


where  P  W  is  the  beam  power  and  r  cm  is  the  beam  radius. 

The  source  term  can  then  be  calculated  by  using  Beer’s  Law  and  multiplying  by  the  absorption  coefficient: 


A  =  naI0e  ^ 


(11-2) 
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Source  Term  (A) 
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Figure  11.1:  Model  calculations  for  a  flat-top  source  term. 


Souce  Term  versus  Axial  Position 


Figure  11.2:  Axial  distribution  of  a  flat-top  source  term. 
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Souce  Term  versus  Radial  Position 


Radial  Position  (cm) 

Figure  1 1.3:  Radial  distribution  of  a  flat-top  source  term. 


At  the  front  of  the  sample  (z  =  0  cm),  A  is  found  to  be  16934.1  W/cm3  and  the  model  predicts  the  same.  Further  into 
the  sample  at  z  —  3  cm,  A  is  found  to  be  313.276  W/cm3  and  the  model  predicts  the  same. 


Source  Term  versus  Radial  Position 


Figure  1 1 .4: 
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Percent  Error  versus  Radial  Position 
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Figure  11.5:  Percent  error  between  model  and  analytical  comparisons  in  figure  11.4 

Figure  11.4  shows  an  analytical  solution  using  Beer’s  Law  to  calculate  the  source  term  versus  the  model  prediction. 
The  largest  percent  error  between  these  curves  is  on  the  order  of  1 0  4  %. 
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11.3  Thermal  Diffusion  Verification 

To  verify  the  thermal  diffusion  in  the  BTEC  model,  the  one-dimensional  code  is  checked  against  known  analytical 
solutions.  Two  analytical  comparisons  are  shown  below. 


11.3.1  One-Dimensional  Stretched  Grid 

The  analytical  solution  in  equation  1 1.3  is  from  [6],  It  is  derived  for  an  infinite  space  (-00  <  x  <  oo).  A  constant  initial 
temperature  (To  =  100  °C)  is  given  for  the  space  (0  <  x  <  L ).  The  analytical  solution  was  calculated  with  L  =  1.0 
cm.  The  BTEC  model  was  run  with  the  constant  temperature  block  from  2.0  cm  to  3.0  cm.  This  is  a  2.5  cm  coordinate 
shift  from  the  analytical  solution.  The  BTEC  model  utilizes  a  stretched  coordinate  grid  to  simulate  an  infinite  problem 
space.  The  coordinate  shift  is  necessary  to  center  the  constant  temperature  block  of  Tq  °C  in  the  center  of  the  uniform 
grid.  The  analytical  solution  is 


T(x,  t ) 
To 


erf[I^=£\  +  erf[L 


\  V4  at 


V  V4at 


where  again 
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(11.3) 


(11.4) 
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Infinite  Space  Validation 


Figure  11.6:  Model  verifications  for  simulating  an  infinite  problem  space 
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Figure  11.7:  Percent  error  for  figure  11.6 


11.3.2  One-Dimensional  Uniform  Grid 

In  this  case,  the  grid  is  uniform  and  extends  from  0.0  cm  to  5.0  cm.  The  entire  problem  space  has  an  initial  temperature 
of  100  °C  and  the  ambient  temperature  is  0  °C.  The  left  boundary  is  an  insulating  boundary  and  the  right  boundary 
is  convective.  The  convective  heat  transfer  rate  is  he  =  .2  W/cm2oC.  The  model  is  compared  to  an  analytical  solution 
form  [6]  after  5  s.  The  analytical  solution  used  here  is  shown  in  equation  1 1.5: 


T(x,t)  = 

m=  1 


Hi  cos(J3mx ) 

L(/32m  +  Hi)  +  Hi  cos(J3mL)  ’ 


(11.5) 
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.8:  Model  verifications  for  convective  and  insulating  boundaries 
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Figure  11.9:  Percent  error  between  model  and  analytical  comparisons  in  figure  11.8 


where  Tq  °C  is  the  initial  temperature  of  the  block  of  size  L  cm  and  Hi  cm'1  is  a  convection  term  and  is  defined  by 


The  f5m  terms  are  the  positive  roots  of 
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/3mtan(fimL)  —  Hi 

and  a  term  is  thermal  diffusivity  and  can  be  found  by 
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(11.6) 


(11.7) 
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The  analytical  solution  1 1.5  used  for  comparison  with  the  BTEC  model  utilized  the  first  100  positive  roots  of  equation 
11.6 


11.3.3  Two-Dimensional  Large  Beam  Verification 

The  one -dimensional  code  simulates  an  infinitely  large  laser  beam.  In  two-dimensional  cases  with  large  beam  diame¬ 
ters,  the  axial  temperature  distribution  of  the  two-dimensional  case  should  compare  to  the  one-dimensional  temperature 
distribution. 
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Figure  11.10:  Temperature  comparison  of  1-D  and  2-D  along  z  axis 
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Figure  11.11:  Percent  error  for  figure  11.10 
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Figures  11.10  and  11.11  show  how  the  axial  temperatures  for  four  different  beam  diameters  in  the  two-dimensional 
model  compare  with  the  one-dimensional  model.  Figure  11.10  shows  the  temperature  distribution  after  3  seconds  and 
figure  11.11  shows  the  temperature  distribution  after  10  seconds.  The  power  density  is  kept  constant  in  all  runs  to 
I  =  3.537  W/cm2.  The  absorption  coefficient  is  pa  =  100  cm-1. 

11.3.4  Two-Dimensional  Mainster  Source 

The  Mainster  source  found  in  [3]  is  a  constant  cylindrical  heat  source  embedded  in  the  problem  space  as  shown  in 
figure  11.12.  The  BTEC  model  has  a  Mainster  source  option.  When  this  option  is  enabled,  the  power,  P  W,  defined 
in  the  emitter  configuration  file  is  used  point  as  the  source  A  W/cm3.  The  cylinder  source  uses  the  beam  diameter 
defined  in  the  emitter  configuration  file  and  is  only  added  to  the  layer  defined  with  its  layer  type  equal  to  1 .  To  produce 
the  embedded  source,  three  layers  are  defined  and  the  center  one  is  defined  as  layer  type  1  from  within  the  layer 
configuration  file. 


Figure  11.12:  Embedded  Mainster  source 


The  analytical  solution  from  11.13  for  this  problem  geometry  is  shown  in  equation  11.8: 
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The  parameters  used  in  the  model  for  comparison  to  the  analytical  solution  in  1 1.8  are  as  follows: 
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All  boundary  conditions  are  considered  sinks  in  the  Mainster  case.  The  temperature  recorded  is  the  maximum  temper¬ 
ature  rise  and  is  found  at  the  center  of  the  heat  source. 
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12.1 

0.3678 

10.0 

14.15 

14.1 

0.3674 

20.0 

15.68 

15.7 

0.1076 

50.0 

17.09 

17.1 

0.0152 

100.0 

17.82 

17.8 

0.1348 

200.0 

18.33 

18.3 

0.1825 

500.0 

18.69 

18.8 

0.5383 

1000.0 

18.77 

19.0 

1.1805 

Table  11.1:  Table  of  values  for  figures  11.13  and  11.14. 


Temperature  Rise  versus  Time 
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Figure  11.13:  Maximum  temperature  rise  comparison  of  analytical  solution  and  BTEC  predictions  at  different  times 
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Percent  Error  versus  Time 


Time  (s) 

Figure  1 1. 14:  Percent  error  of  BTEC  as  compared  to  analytical  solution  with  Mainster  source  term 


11.4  COMSOL  Thermal  Diffusion  Comparison 


The  BTEC  model  data  was  compared  to  COMSOL  data  provided  by  Dr.  Irwin  Goldberg  and  Misty  Garcia.  The 
COMSOL  model  was  a  two-dimensional  Cartesian  model  with  a  circular  source  in  the  center.  The  source  was  a  5  mm 
beam  that  provides  100  W/cm3  as  a  source  term  for  the  model.  A  constant  source  in  a  two-dimensional  model  is  the 
same  as  and  infinite  source  in  a  three-dimensional  model.  The  BTEC  was  configured  as  a  constant  source  from  z  =  0 
cm  to  z  =  10  cm.  Only  points  in  the  z  =  5  cm  plane  were  considered.  The  maximum  radial  boundary  was  simulated 
to  be  infinity  as  well.  Several  points  were  compared  between  the  two  models.  Figure  11.15  shows  two  points  and  how 
they  compare  between  the  two  models. 


0123456789  10 

Time  [s] 

Figure  11.15:  Infinite  source  COMSOL  comparison 


The  models  were  compared  in  figure  1 1. 15  at  points  2.4  mm  and  3.4  mm  from  the  center  of  the  source.  The  first  point 
would  be  slightly  inside  the  source  and  the  second  is  outside  the  source. 
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11.5  Damage  Integral 


To  validate  the  BTEC  implementation  of  the  Arrhenius  damage  integral  (equation  9.1),  a  separate  model  was  created 
in  Matlab.  The  model  used  an  array  of  temperature  versus  time  values  to  numerically  compute  the  value  of  the  damage 
integral  at  each  time  included  in  the  array.  Temperature  versus  time  and  damage  versus  time  arrays  were  created  for  a 
specific  grid  point  using  the  BTEC  model.  The  temperature  versus  time  array  was  used  in  the  Matlab  model  to  create 
a  damage  versus  time  array  which  is  compared  to  the  BTEC’s  damage  versus  time  array  in  figure  11.16. 


Damage  versus  Time 


Figure  11.16:  Probability  of  damage  comparisons  between  BTEC  and  Matlab  models 


This  comparison  was  done  with  a  constant  laser  source  of  50  W  with  a  beam  diameter  of  0.1  cm.  The  material  had  the 
following  parameters: 
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11.6  Beam  Propagation 


To  help  verify  that  the  beam  propagation  methods  in  the  BTEC  model  are  working  correctly,  three  separate  methods 
are  compared  for  the  same  case.  A  Gaussian  beam  is  focused  for  1  cm  using  the  GBP,  Hankel  transform,  and  Pade 
beam  propagation  methods  and  the  results  are  compared  in  figure  11.17. 
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Irradiance  versus  Radial  Position 
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Figure  11.17:  Focused  Gaussian  comparison  between  Hankel,  Pade,  and  GBP 

To  accomplish  the  comparison  in  figure  1 1.17,  the  GBP  model  is  configured  with  a  spherical  wave  front  at  z  =  0  to  be 
propagated  by  either  the  Hankel  or  Pade  method  into  the  sample.  After  initial  data  are  generated  (before  any  heating) 
with  the  Pade  and  Hankel  method,  data  are  extrapolated  for  irradiance  at  all  radial  positions  for  z  =  1  cm.  The  model 
is  then  configured  to  have  the  initialize  GBP  a  plane  at  z  =  1  cm.  The  GBP  computes  the  translation  of  the  Gaussian 
beam  to  the  z  =  1  cm  plane. 
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Chapter  12 


Configuration  and  Use 


12.1  Overview 

The  BTEC  thermal  model  is  a  2-D  (cylindrical  coordinate  system)  simulation  of  optical-tissue  thermal  interactions. 
The  code  supports  the  illumination  of  tissues  by  optical  sources,  the  temperature  response  from  the  linear  absorption 
of  optical  radiation,  and  subsequent  damage. 

An  alternating  direction  implicit  (ADI)  finite  difference  method  is  employed  in  the  solution  of  the  heat  equation. 
The  method  is  commonly  referred  to  as  the  Peaceman-Rachford  method.  The  simulation  can  be  configured  with  a 
source  term  defined  by  a  single  or  multiple  optical  emitters  (laser  or  broadband).  The  tissue  simulation  is  represented 
as  a  one-dimensional  stack  of  homogeneous  layers  along  the  z-axis.  Each  layer  can  have  differing  thermal,  optical,  and 
physical  properties.  The  linear  absorption  coefficient  of  each  layer  defines  the  energy  transfer  from  the  optical  source 
to  the  tissue.  Boundary  conditions  along  the  axial  direction  may  be  selected  as  a  sink  (temperature  held  constant)  or 
as  a  combination  of  surface  boundary  conditions  (convection,  radiative,  or  evaporative). 

The  model  is  configured  through  three  types  of  configuration  files:  (1)  The  main,  or  top-level  configuration  file  (con- 
fig.*. in)  which  specifies  the  problem  geometry,  simulation  time,  boundary  condition  types,  and  search  parameters 
used  to  find  tissue  damage  thresholds.  (2)  Emitter  files  (*. emitter)  which  specify  the  optical  emitter  properties  such  as 
wavelength  or  spectrum,  power,  spatial  profile,  time-dependent  power,  and  parameters  which  specify  time-step  size  in 
the  simulation.  (3)  Layer  files  (*. layer)  which  specify  thermal,  optical,  and  physical  properties  of  homogeneous  tissue 
layers  along  with  the  thickness  of  each  layer. 

Running  BTECthermal  with  the  options  -c  will  generate  template  versions  of  all  the  configuration  files.  The  tem¬ 
plates  are  listed  here. 


12.2  Top-Level  (config.*.in)  Configuration 

This  is  the  top-level  configuration  file  for  the  BTEC  Thermal  Model.  This  is  the  file  where  the  user  inputs  the  desired 
overall  simulation  parameters: 

The  following  is  a  short  description  of  each  input: 

1 .  #Key  Value  header:  The  first  line  of  the  config  must  start  with  the  string  #Key  Value. 

2.  Dimensions:  (INTEGER,  1  \2)  The  number  of  dimensions  to  use  in  the  simulation.  Currently,  BTEC  supports 
1-D,  and  2-D  cylindrical  coordinates. 
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3.  SimulationType:  (INTEGER,  0  |2)  This  option  is  depreciated  and  should  be  set  to  0.  Previous  versions  of 
BTEC  separated  Propagation  emitters  from  the  others.  Propagation  emitters  are  now  implemented  through  the 
Standard  Emitter  interface,  and  specified  by  the  StandardEmitter  entries  below.  However,  z-scan  simulations 
are  still  invoked  from  here. 

4.  AxialGridType:  (INTEGER,  0  1 1 )  The  user  can  choose  to  use  a  uniform  grid  (0),  or  stretched  grid  ( 1 )  in  the  axial 
direction.  In  a  uniform  grid,  the  spacing  between  values  on  the  axis  is  constant.  The  stretched  grid  adds  points 
to  the  end(s)  of  the  uniform  grid,  and  spacing  between  these  extra  points  grows.  Grid  stretching  only  occurs  at 
the  end  of  the  axis,  and  will  only  be  used  when  the  boundary  condition  is  set  to  a  SINK  for  that  end.  It  will  not 
be  used  if  the  end  of  the  axis  represents  a  surface  boundary.  An  example  of  a  stretched  grid  used  by  Mainster  et 
al.  is  given  on  pg  403  of  Welch  and  Gemert’s  book,  Optical-Thermal  Response  of  Laser-Irradiated  Tissue. 

5.  Nz:  (INTEGER,  >1)  The  number  of  axial  divisions.  The  axial  grid  will  consist  of  this  number  plus  one  grid 
points. 

6.  zMin:  (DOUBLE  PRECISION,  >0.0)  The  minimum  axial  coordinate  value.  Units  are  centimeters  (cm).  This 
coordinate  applies  to  the  uniform  portion  of  the  grid. 

7.  zMax:  (DOUBLE  PRECISION,  >0.0)  The  maximum  axial  coordinate  value.  Units  are  centimeters  (cm).  This 
coordinate  applies  to  the  uniform  portion  of  the  grid. 

8.  zStretchRatio:  (DOUBLE  PRECISION,  >1.0)  The  grid  spacing  stretching  ratio  for  the  extended,  stretched  axial 
grid.  Each  A z  calculated  by  this  ratio  times  the  previous  A’.  The  value  must  be  relatively  small  for  numerical 
derivative  approximations  to  hold,  and  to  limit  the  extent  of  the  grid.  This  parameter  is  ignored  for  axial  grid 
type  of  0  (uniform  grid),  and  when  the  minimum  or  maximum  axial  coordinate  boundary  condition  is  set  to  a 
surface  boundary  condition  type. 

9.  zMinBC:  (INTEGER,  0-7)  The  axial  grid  may  be  assigned  either  a  sink  boundary  condition  mechanisms).  Cur¬ 
rently,  the  model  supports  the  following  assignments  for  the  Minimum  (and  Maximum)  Axial  Position  Boundary 
Condition.  Specific  parameters  for  the  boundary  conditions  are  specified  within  the  layer  configuration  files. 


TYPE  VALUE 

SINK  0 

CONVECTIVE  1 

RADIATIVE  2 

CONVECTIVE  +  RADIATIVE  3 

EVAPORATIVE  4 

CONVECTIVE  +  EVAPORATIVE  5 

RADIATIVE  +  EVAPORATIVE  6 

RADIATIVE  +  EVAPORATIVE  +  CONVECTIVE  7 


10.  zMaxBC:  (INTEGER,  0-7)  See  parameter  9  for  a  description,  the  same  condition  types  apply. 

11.  RadialGridType:  (INTEGER,  0  1 1)  If  the  user  selects  2  dimensions,  then  the  second  axis  (r)  can  be  stretched 
at  the  maximum  end.  Surface  boundary  conditions  are  not  implemented  at  the  rmax  boundary,  so  stretching  can 
always  be  used  when  specified  here. 

12.  Nr:  (INTEGER,  >1)  The  number  of  radial  divisions.  The  radial  grid  will  consist  of  this  number  plus  one  grid 
points. 

13.  rMax:  (DOUBLE  PRECISION,  >0.0)  The  maximum  radial  coordinate,  units  are  centimeters  (cm).  This  coordi¬ 
nate  applies  to  the  uniform  portion  of  the  grid.  If  a  stretched  grid  is  selected  then  additional  grid  points  are  added 
to  the  grid  with  incremental  spacing  ratios  as  set  below.  (Note:  The  minimum  radial  coordinate  is  assumed  to 
be  r  —  0.) 


78 

DISTRIBUTION  A.  APPROVED  FOR  PUBLIC  RELEASE 


14. 


15. 

16. 

17. 


18. 


19. 

20. 


21. 


22. 

23. 

24. 


25. 


26. 


27. 


rStretchRatio:  (DOUBLE  PRECISION,  >  1 .0)  The  grid  spacing  stretching  ratio  for  the  extended,  stretched  radial 
grid.  Each  grid  point  is  spaced  by  this  ratio  times  the  spacing  of  the  previous  grid  point  spacing.  The  value  must 
be  relatively  small  for  numerical  derivative  approximations  to  hold,  and  to  limit  the  radial  extent  of  the  grid. 
This  parameter  is  ignored  for  radial  grid  type  of  0  (uniform  grid). 


rMaxBC:  (INTEGER,  0  |1)  The  radial  grid  may  be  assigned  either  a  sink  boundary  condition  (A T(rmax)  =  0)  or 
an  insulating  boundary  condition  (  %\r_r  =  o). 


TotalSimTime:  (DOUBLE  PRECISION,  >0.0)  The  simulation  time  starts  at  t  -  0.0  and  runs  to  this  user- 
specified  value,  given  in  seconds. 


dt:  (DOUBLE  PRECISION,  >0.0)  The  simulation  will  make  time  steps  of  this  value,  unless  an  adaptive  time 
step  is  specified  in  the  first  emitter  configuration  file.  In  the  case  of  an  adaptive  time  step,  this  parameter  is 
ignored.  Ideally,  the  time  step  should  be  set  to  a  value  much  shorter  than  the  thermal  relaxation  time,  such  that 
the  numerical  estimates  for  derivatives  remain  accurate. 


dtMax:  (DOUBLE  PRECISION,  >0.0)  In  the  case  of  adaptive  time  steps,  the  maximum  time  step  is  currently 
not  bounded  by  an  estimate  for  stability.  This  parameter  allows  the  user  to  limit  the  time  step  to  a  maximum 
value  if  conditions  of  instability  or  inaccuracy  are  observed  the  solution. 

TissueBaselineTemp:  (DOUBLE  PRECISION,  <0.0)  This  is  the  temperature,  in  °C,  of  the  tissue  corresponding 
to  AT  equal  to  zero  in  the  heat  equation  solution.  The  value  is  typically  37.0  for  body  temperature. 

AmbientTemp:  (DOUBLE  PRECISION,  <0.0)  This  value  is  the  temperature  of  the  surrounding  atmosphere  for 
conditions  where  surface  boundary  conditions  are  employed.  The  value  must  be  present  but  will  not  be  used  if 
both  the  zmin  and  zmax  boundary  conditions  are  set  to  SINK. 

RelHumidity:  (DOUBLE  PRECISION,  0.0  -  1 .0)  This  value  is  the  relative  humidity  value  for  the  atmosphere 
surrounding  the  sample  and  is  employed  only  for  surface  boundary  conditions  using  the  evaporative  boundary 
condition  type,  or  combinations  of  surface  boundary  conditions  using  the  evaporative  boundary  condition  type. 

LogDataFlag:  (INTEGER,  0  |1)  This  turns  logging  on  and  off.  Several  hies  are  generated  when  this  Hag  is  set 
to  1,  such  as  the  thermal  distribution  at  some  time,  the  source  term  at  some  time,  and  so  on. 

Loglnterval:  (INTEGER,  >0)  This  integer  specifies  the  number  of  time  steps  between  logging.  Note  that  in  the 
case  of  the  non-uniform  time  step  the  logged  data  will  not  be  equally  spaced  in  time. 

UserDamageThreshold:  (DOUBLE  PRECISION,  >0.0)  The  valued  used  in  the  evaluation  of  the  damage  inte¬ 
gral  to  indicated  that  tissue  has  been  damaged.  The  damage  integral  is  computed  at  each  grid  point  during  the 
simulation.  If  the  damage  value  becomes  greater  than  this  value,  then  the  tissue  is  flagged  as  damaged. 

DamageThresholdSearchFlag:  (INTEGER,  0  |1)  This  sets  the  option  for  the  simulation  to  search  for  a  threshold 
power  to  give  a  user-specified  damage  condition.  If  the  user  would  rather  find  the  damage  given  by  a  certain  set 
power  and  not  do  a  search,  this  would  be  set  to  0.  Entering  1  indicates  that  a  search  should  be  done. 

MinPowerRatio/MaxPowerRatio:  (DOUBLE  PRECISION,  >0.0)  These  inputs  are  the  multipliers  of  the  power 
specified  for  the  emitter.  They  give  the  simulation  bounds  with  which  to  start  the  search.  If  the  threshold  power 
is  outside  of  this  range,  the  simulation  will  print  an  error  message  and  exit. 

ThreshSearchType:  (INTEGER,  0-3)  There  are  four  different  search  types.  A  search  type  of  0  means  search 
for  the  threshold  to  cause  damage  at  any  point  within  the  computational  grid.  A  search  type  of  1  means  to  look 
for  a  specific  temperature  rise  at  any  point  within  the  grid.  A  search  type  of  2  means  to  look  for  damage  of  a 
given  radius.  A  search  type  of  3  means  to  look  for  a  given  temperature  rise  out  to  a  given  radius.  For  cases 
where  damage  is  being  examined,  the  criterion  for  damage  is  that  the  damage  integral  is  greater  than  or  equal  to 
1. 
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28.  PrimeThresh Value:  (DOUBLE  PRECISION,  <0.0)  This  value  depends  on  the  search  type.  For  a  search  type  of 
0,  this  parameter  is  not  used.  For  a  search  type  of  1  or  3,  this  is  the  temperature  rise.  For  a  search  type  of  2,  this 
is  the  radius  of  damage  the  must  be  present  as  a  threshold  criteria.  (Note:  At  the  current  model  ability  damage, 
area  should  be  equal  or  greater  to  the  beam  area  to  prevent  variations  in  thresholds  with  a  change  in  the  number 
of  grid  points.) 

29.  SecThreshValue:  (DOUBLE  PRECISION,  >0.0)  This  number  is  only  used  for  search  type  3.  It  is  the  radius  at 
which  the  temperature  rise  specified  must  occur. 

30.  ConvergThresh:  (DOUBLE  PRECISION,  >0.0)  This  is  the  percent  error  in  threshold  criterion  that  the  user 
will  allow  in  their  results.  The  code  will  iterate  until  the  predicted  threshold  criterion  is  within  this  ratio  of 
agreement,  and  stops  the  search. 

31.  Emitterfi]:  (STRING,  i  =  0, 1,2, ...)  The  string  specifies  the  path  name  (relative  or  full)  to  the  Emitter  files.  The 
index  i  must  start  at  0,  and  increases  for  each  Emitter. 

32.  StandardEmitterfi]:  (STRING,  i  —  0, 1,2, ...)  The  string  specifies  the  path  name  (relative  or  full)  to  the  Standard- 
Emitter  files.  Propagation  emitters  and  Monte  Carlo  Scattering  emitters  are  implemented  through  this  interface. 
The  index  i  must  start  at  0,  and  increases  for  each  Emitter. 

33.  Layer[/]:  (STRING,  i  =  0, 1,2, ...)  The  string  specifies  the  path  name  (relative  or  full)  to  the  Layer  files.  The 
index  i  must  start  at  0,  and  increases  for  each  Layer. 

34.  InitialConditionsFlag:  (INTEGER,  0  1 1 )  This  flag  indicates  that  the  user  has  specified  the  loading  of  an  initial 
temperature  distribution  within  the  solution  grid.  The  initial  condition  is  specified  by  the  file  below. 

35.  InitialConditionsFile:  (STRING)  The  file  containing  the  initial  conditions  temperature  profile. 

36.  PropagationType:  (INT)  This  configuration  is  obsolete.  See  StandardEmitter  configuration  for  Propagation 
emitter  configuration. 

37.  NumPropSteps:  (INT)  This  configuration  is  obsolete.  See  StandardEmitter  configuration  for  Propagation  emitter 
configuration. 

38.  ApertureCoordinate:  (DOUBLE  PRECISION)  This  specifies  the  z  coordinate  of  the  aperture  for  z-scan  simula¬ 
tions. 

39.  ApertureRadius:  (DOUBLE  PRECISION,  >0)  The  aperture  radius  for  z-scan  simulations. 

40.  NDiffPattern:  (INT,  >0)  The  number  of  points  to  sample  in  the  diffraction  pattern  calculations. 

41.  DiffPatternMaxRadius:  (DOUBLE  PRECISION,  >0.0)  Maximum  radial  coordinate  in  diffraction  pattern  calcu¬ 
lations. 

42.  ZScanStartCoord:  (DOUBLE  PRECISION,  <ApertureCoordinate)  The  initial  sample  position  for  z-scan  simu¬ 
lation 

43.  ZScanEndCoord:  (DOUBLE  PRECISION,  < ApertureCoordinate)  The  final  sample  position  for  z-scan  simula¬ 
tion 

44.  ZScanStepSize:  (DOUBLE  PRECISION,  <ApertureCoordinate)  Step  size  taken  between  initial  and  final  sam¬ 
ple  positions  in  z-scan  experiment. 

45.  ZScanLensCoord:  (DOUBLE  PRECISION,  <ZScanStartCoord)  Lens  coordinate  for  z-scan  simulations.  Lens 
must  be  placed  in  front  of  all  sample  positions. 

46.  ZScanLensFocalLength:  (DOUBLE  PRECISION)  The  focal  length  of  the  lens  for  z-scan  simulations. 
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12.3  Emitter  (Remitter)  Configuration 


1.  EmitterType:  (INTEGER,  1  |2)  Set  to  1  for  laser  source,  2  for  broadband  source. 

2.  PulseType:  (INTEGER,  1-3)  This  determines  of  the  emitter  produces  a  single  (1),  square  pulse  ;  multiple  (2) 
square  pulses;  or  a  user  defined  (3)  pulse  train. 

3.  FocusType:  (INTEGER,  1  |2)  Specifies  if  the  focus  is  stationary  (1)  at  some  z  position,  or  moves  (2). 

4.  TimeStepType:  (INTEGER,  1  |2)  The  time  step  can  either  be  fixed  (1),  with  a  constant  dt,  or  adaptive  (2),  where 
dt  grows  for  each  time  step1 

5.  ProfileType:  (INTEGER,  1  |2)  The  beam  profile  can  be  set  to  a  Gaussian  (1),  Flattop  (2),  Annular  (3),  or  User 
defined  (4). 

6.  PeakPower:  (DOUBLE  PRECISION,  >=  0.0)  The  peak  power  of  the  emitter  in  W. 

7.  MinWavelength:  (DOUBLE  PRECISION,  >=  0.0)  For  a  laser  source,  this  is  the  wavelength  of  the  laser.  For  a 
broadband  source,  this  specifies  the  minimum  wavelength  in  the  source. 

8.  MaxWavelength:  (DOUBLE  PRECISION,  >=  MinWavelength)  For  a  Broad  Band  source,  this  is  the  maximum 
wavelength  present  in  the  source. 

9.  DeltaWavelength:  (DOUBLE  PRECISION,  >0.0)  For  a  Broad  Band  source,  the  source  is  assumed  to  contain  a 
set  of  wavelength  between  MinWAvelength  and  MaxWavelength,  with  each  wavelength  having  a  power  associ¬ 
ated  with  it.  This  sets  the  resolution  for  the  sampling  of  the  wavelengths  in  the  source. 

10.  BeamDiameter:  (DOUBLE  PRECISION,  >0.0)  The  beam  diameter.  For  a  Gaussian  beam,  this  specifies  the 
1/e2  diameter. 

11.  BeamWaistPosition:  (DOUBLE  PRECISION,  >0.0)  Z-coordinate  of  emitters  focus.  Beam  will  be  focused  at 
this  coordinate. 

12.  BeamDivergence:  (DOUBLE  PRECISION)  The  divergence  of  the  beam,  at  the  focus. 

13.  PulseDuration:  (DOUBLE  PRECISION,  >0.0)  The  temporal  pulse  width. 

14.  PulsePeriod:  (DOUBLE  PRECISION,  >PulseDuration)  The  time  between  the  front  of  multiple  pulses. 

15.  StartTime:  (DOUBLE  PRECISION,  >=  0.0)  The  time  that  the  emitter  comes  on  and  starts  pulsing. 

16.  StopTime:  (DOUBLE  PRECISION,  >StartTime)  The  time  that  the  emitter  turns  off. 

17.  dtMinOn:  (DOUBLE  PRECISION,  >0.0)  The  minimum  time  step  that  the  adaptive  time  step  will  take  when 
the  emitter  it  turned  on.  This  is  what  the  time  step  is  set  to  when  an  emitter  is  switched  on. 

18.  dtMinOff:  (DOUBLE  PRECISION,  >0.0)  The  minimum  time  step  that  the  adaptive  time  step  will  take  when 
the  emitter  it  turned  off.  This  is  what  the  time  step  is  set  to  when  an  emitter  is  switched  off. 

19.  StretchOn:  (DOUBLE  PRECISION,  >1.0)  The  stretch  ratio  used  for  the  adaptive  time  step  when  the  emitter  is 
on. 

20.  StretchOff:  (DOUBLE  PRECISION,  >1.0)  The  stretch  ratio  used  for  the  adaptive  time  step  when  the  emitter  is 
off. 

ldt  starts  at  a  minimum  time  step  and  grows  with  each  time  step  until  an  emitter  state  changes.  So  if  an  emitter  is  off  at  some  time,  and  turns  on 
during  the  next  time  step,  dt  is  reset  to  the  minimum  value. 


81 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


12.4  Standard  Emitter  (Remitter)  Configuration 


Standard  emitters  are  indicated  in  the  main  configuration  file  by  the  StandardEmitter[i]  tag.  Many  of  the  same  con¬ 
figuration  options  listed  in  the  emitter  configuration  also  apply  for  the  standard  emitter  configuration.  The  parameters 
specific  to  the  standard  emitters  are  listed  below. 


1.  EmitterType:  (STRING)  The  type  in  the  standard  emitter  configuration  is  given  by  a  string.  Currently  supported 
strings  are  ”Pade”,  ’’Hankel”,  ’’Scattering”. 

2.  PhotonPacketDensity:  (DOUBLE  PRECISION,  >0.0)  For  ’’Scattering”  EmitterType.  The  number  of  photons  to 
launch  in  a  1  cm2  area. 

3.  Seed:  (INTEGER,  >0)  For  ’’Scattering”  EmitterType.  Seed  for  initializing  random  number  generator. 

4.  NoiseReductionFactor:  (INTEGER,  >=  0)  For  ’’Scattering”  EmitterType.  Used  for  smoothing  the  source  term 
for  a  scattering  emitter  near  the  r  -  0  axis. 

5.  zResFactor:  (INTEGER,  >=  1)  For  ”Pade”  EmitterType.  Resolution  factor  for  Pade  z  axis.  Pade  propagation 
grid  requires  higher  resolution  than  the  thermal  grid.  The  propagation  grid  resolution  is  based  on  the  thermal 
grid  resolution,  multiplied  by  a  resolution  factor. 

6.  rResFactor:  (INTEGER,  >=  1)  For  ’’Pade”  EmitterType.  Resolution  factor  for  Pade  r  axis.  Pade  propagation 
grid  requires  higher  resolution  than  the  thermal  grid.  The  propagation  grid  resolution  is  based  on  the  thermal 
grid  resolution,  multiplied  by  a  resolution  factor. 

7.  OpticsFile:  (STRING)  For  ’’Pade”  and  ’’Hankel”  EmitterTypes.  If  given,  defines  an  optics  setup  for  the  Gaussian 
Beam  Propagation.  Currently,  this  should  be  set  to  ’’EYE”,  which  gives  an  initial  wavefront  for  the  propagation 
emitter  corresponding  to  the  curvature  of  the  eye. 

8.  P:  (INTEGER,  >=  0)  For  ’’Hankel”  EmitterType.  The  order  of  Bessel  Functions  to  use  in  Hankel  propagation. 
This  should  be  set  to  zero. 

9.  numBesselFunctions:  (INTEGER,  0  -  4999)  For  ’’Hankel”  EmitterType.  The  number  of  Bessel  Frequencies  to 
use  in  Hankel  transforms. 

10.  BesselPrecision:  (STING,  ’’single”  |”double”)  For  ’’Hankel”  EmitterType.  Precision  to  use  in  Bessel  Function 
evaluation.  Double  precision  evaluation  can  require  much  more  computational  time. 

12.5  Layer  (*.layer)  Configuration 

1 .  LayerType:  (INTEGER)  An  integer  that  tags  the  layer. 

2.  Description:  (STRING)  A  description  of  the  layer. 

3.  Thickness:  (DOUBLE  PRECISION,  >0.0)  The  layer  thickness 

4.  Density:  (DOUBLE  PRECISION,  >0.0)  The  layer  density,  given  in  g/cm 3. 

5.  SpecificHeat:  (DOUBLE  PRECISION,  >0.0)  Specific  heat  of  the  material. 

6.  Conductivity:  (DOUBLE  PRECISION,  >0.0)  Conductivity  of  the  material. 

7.  ConvHeatTransRate:  (DOUBLE  PRECISION,  >0.0)  The  convective  heat  loss  from  an  air  boundary  with  the 
material. 

8.  Emissivity:  (DOUBLE  PRECISION,  0.0  -  1.0)  The  emissivity  of  the  material. 
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9.  BloodFlowRate:  (DOUBLE  PRECISION,  >=  0.0)  The  blood  flow  rate,  given  in  glow’s,  through  the  material. 
Setting  this  to  0.0  turns  off  blood  flow. 

10.  Refractivelndex[/]:  (STRING)  A  string  specifying  the  refractive  index  (n)  and  change  in  refractive  index  (4p0 
for  a  given  wavelength.  Values  for  multiple  wavelength  can  be  given  by  indexing  with  i.  in  the  STRING, 
wavelength  is  given  first,  then  index  of  refraction,  then  change  in  index  of  refraction  with  temperature.  So,  for 
example,  if  the  index  of  refraction  for  the  material  at  2000  nm  is  1.35  and  it  changes  with  temperature  at  a 
constant  rate  of  0. 1,  the  entry  would  be  given  as: 

RefractiveIndex[0]  =  ”2000  1.35  0.1” 

11.  Absorption[/]:  (STRING)  A  string  specifying  the  absorption  coefficient  of  the  material  for  a  given  wavelength. 
Again,  multiple  wavelengths  can  be  given  yy  indexing  i. 

12.  Anisotropy [/]:  (STRING)  A  string  specifying  the  anisotropy  of  the  material  at  a  wavelength.  Multiple  wave¬ 
lengths  can  be  given. 

13.  Scattering [i]:  (STRING)  A  string  specifying  the  scattering  coefficient  of  the  material  at  a  wavelength.  Multiple 
wavelengths  can  be  given. 

14.  Reflectance  [/]:  (STRING)  A  string  specifying  the  reflectance  of  the  material  at  a  wavelength.  Multiple  wave¬ 
lengths  can  be  given. 

15.  Temp[i],  A[z],  Ea[z]  (DOUBLE  PRECISION)  These  values  specify  the  normalization  constant  (A)  and  the  rate 
process  coefficient  ( Ea )  for  the  damage  integral  in  the  material.  When  calculating  the  damage  integral,  the  sim¬ 
ulation  will  look  for  the  two  temperatures  (given  in  Kelvin)  that  bracket  the  current  temperature  and  use  the 
coefficients  given  for  the  lower  bound.  For  example,  if  the  following  coefficients  are  listed  in  the  layer  configu- 


ration  file: 

Temp[0] 

=  250 

A[0] 

=  10E45 

Ea[0] 

=  6.28E10 

Temp[l] 

=  300 

A[l] 

=  10E99 

Ea[l] 

=  3.33E20 

Temp  [2] 

=  500 

A[2] 

=  10E130 

Ea[2] 

=  2.08E23 

and  the  current  temperature  is  350  K,  then  the  two  coefficients  that  will  be  used  are  A  =  10"  and  Ea  =  3.3320 
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#KeyValue  config  file 

Dimensions 

SimulationType 

AxialGridType 

Nz 

zMin 

zMax 

zStretchRatio 

zMinBC 


=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 


zMaxBC 


=  VALUE 


RadialGridType 

Nr 

rMax 

rStretchRatio 

rMaxBC 


=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 


# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 


INT: 
INT: 
INT: 
INT: 
DBL : 
DBL : 
DBL: 
INT: 


INT: 


INT: 

INT: 

DBL: 

DBL: 

INT: 


number  of  dimensions  in  simulation  (1  |  2) 

problem  (sim)  type  (0  STD  THERMAL,  1  PROPAGATION+THERMAL,  2  Z-SCAN) 
0  =  UNIFORM,  1  =  STRETCHED 
number  of  z  grid  divisions 


minimum  z  coordinate 
maximum  z  coordinate 

(**)  z  stretch  ratio  for  non-uniform  grid 

z  min  boundary  condition  type 

SINK  0 

CONVECTIVE  1 

RADIATIVE  2 

CONVECTIVE  +  RADIATIVE  3 

EVAPORATIVE  4 

CONVECTIVE  +  EVAPORATIVE  5 

RADIATIVE  +  EVAPORATIVE  6 

ALL  SURFACE  7 

z  min  boundary  condition  type 

SINK  0 

CONVECTIVE  1 

RADIATIVE  2 

CONVECTIVE  +  RADIATIVE  3 

EVAPORATIVE  4 

CONVECTIVE  +  EVAPORATIVE  5 

RADIATIVE  +  EVAPORATIVE  6 

ALL  SURFACE  7 


radial  grid  type  (0  =  UNIFORM,  1  =  STRETCHED) 
number  of  r  grid  divisions 
Max  r  coordinate 

(**)  r  stretch  ratio  for  non-uniform  grid 
r  max  boundary  condition  type  (0=SINK;  l=INSULATOR) 


TotalSimTime 

dt 

dtMax 

TissueBaseTemp 
Ambient Temp 
Re  1 Humidity 
LogDataFlag 
Loglnterval 


=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 


#  DBL:  total  simulation  time 

#  DBL:  time  step  for  fixed  emitter,  not  for  adaptive 

#  DBL:  maximum  time  step  for  adaptive  emitter  (for  stability) 

#  DBL:  tissue  baseline  temperature 

#  DBL:  ambient  temperature 

#  DBL:  relative  humidity  (0-1.0,  .6  =  60%) 

#  INT:  log  data  flag  (0  =  NO,  1=YES) 

#  INT:  log  interval  in  time  steps  (10  =  LOG  1  in  10) 


UserDamageThresh 

DamageThreshSearchFlag 

MaxPowerRatio 

MinPowerRatio 

ThreshSearchType 

PrimeThreshValue 

SecThreshValue 

ConvergThresh 


VALUE 

VALUE 

VALUE 

VALUE 

VALUE 

VALUE 

VALUE 

VALUE 


#  DBL:  User-defined  damage  threshold  value  (damage  integral  value) 

#  INT:  damage  threshold  search  flag  (0  =  NO,  1  =  YES  ) 

#  DBL:  (##)  max  power  ratio  for  threshold  search 

#  DBL:  (##)  min  power  ratio  for  threshold  search 

#  INT:  (##)  threshold  search  type  (0=damage, l=dT, 2=Les  Rad, 3=rad  dT) 

#  DBL:  (##)  primary  threshold  value  (radius,  dT,  radius  dT) 

#  DBL:  (##)  secondary  threshold  value  (dT  in  Thresh  Search  Type  3) 

#  DBL:  (##)  convergence  threshold  for  search  (e.g.  0.05  =  5%  ) 


Emitter [0] 
Layer [0] 
Layer [1] 


=  "STRING"  #  STR:  emitter  file 
=  "STRING"  #  STR:  layer  filename 
=  "STRING"  #  STR:  layer  filename 


InitialConditionsFlag  =  VALUE  #  INT:  initial  conditions  flag  (0  =  none,  1  =  read  initial  conditions  file) 

InitialConditionsFile  =  "STRING"  #  STR:  initial  conditions  file 


PropagationMethod 

=  VALUE 

# 

INT: 

propagation  method  (0=FDBPM,  1=PADE,  2=PADE,  Non-Uniform  grid) 

NumPropSteps 

=  VALUE 

# 

INT: 

(++) 

number  of  prop  steps  per  grid  interval 

AperatureCoordinate 

=  VALUE 

# 

DBL: 

(++) 

aperture  or  screen  coordinate  [cm] 

AperatureRadius 

=  VALUE 

# 

DBL: 

(++) 

aperture  radius  [cm] 

NDif fPattern 

=  VALUE 

# 

INT: 

(++) 

diffraction  pattern  num  points 

Dif fPatternMaxRadius 

=  VALUE 

# 

DBL: 

(++) 

diffraction  pattern  max  radius 

ZScanSt art Coord 

=  VALUE 

# 

DBL: 

(++) 

z-scan  start  coordinate  [cm] 

ZScanStopCoord 

=  VALUE 

# 

DBL: 

(++) 

z-scan  end  coordinate  [cm] 

ZScanStepSize 

=  VALUE 

# 

DBL: 

(++) 

z-scan  coordinate  step  size  [cm] 

ZScanLensCoord 

=  VALUE 

# 

DBL: 

(++) 

z-scan  lens  coordinate  [cm] 

ZScanLensFocalLength 

=  VALUE 

# 

DBL: 

(++) 

z-scan  lens  focal  length  [cm] 

//-  Notes: 

//-  (**)  Not  used  for  uniform  grid  spacing  —  code  sets  to  1.0 

//-  (##)  Only  Used  for  Threshold  Searches 

//-  (++)  Only  Used  for  Propagation  Analysis  (and  secondary  prop  grid) 


//-  dt  is  set  by  emitter  when  variable  timestep  emitters  are  used 

//-  See  documentation  in  /doc  dirctory  for  additional  parameter  information 


Figure  12.1:  Example  of  the  main  configuration  file 
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#KeyValue  emitter  configuration 


EmitterType 

= 

VALUE 

PulseType 

= 

VALUE 

FocusType 

= 

VALUE 

TimeStepType 

= 

VALUE 

SpectrumType 

= 

VALUE 

Prof ileType 

= 

VALUE 

PeakPower 

= 

VALUE 

MinWaveLength 

= 

VALUE 

MaxWaveLength 

= 

VALUE 

DeltaWaveLength 

= 

VALUE 

BeamDiameter 

= 

VALUE 

BeamWaistPosition 

= 

VALUE 

BeamDivergence 

= 

VALUE 

PulseDuration 

= 

VALUE 

PulsePeriod 

= 

VALUE 

StartTime 

= 

VALUE 

StopTime 

= 

VALUE 

dtMinOn 

= 

VALUE 

dtMinOf f 

= 

VALUE 

StretchOn 

= 

VALUE 

StretchOf f 

= 

VALUE 

SARFilename 

= 

"STING" 

PulseFilename 

= 

"STING" 

PowerSpectrumFilename 

= 

"STING" 

Prof ileFilename 

= 

"STING" 

FocusFilename 

= 

"STING" 

#  INT :  emitter  type  (1  Lsr,  2  BB,  3  Mnstr,  4  Prop,  5  SAR,  6  MCML) 

#  INT:  pulse  type  (1  Single,  2  Multiple,  3  USER) 

#  INT:  focus  type  (1  Fixed,  2  Moving  [thermal  lens]) 

#  INT:  time  step  type  (1  Fixed,  2  Adaptive) 

#  INT:  for  broad  band  source,  type  of  wavelength  spectrum  (1  Black  Body,  2  USER) 

#  INT:  profile  type  (1  Gaussian,  2  Tophat,  3  Annular,  4  USER) 

#  DBL:  peak  power  of  emitter  [W] 

#  DBL:  minum  wavelength  for  broad  band  source,  or  wavelength  of  laser  source  [1/cm] 

#  DBL:  minum  wavelength  for  broad  band  source 

#  DBL:  wavelength  step  to  take  for  numerical  integration 

#  DBL:  beam  diameter  [cm]  Note:  this  is  1/e  for  gaussian  beam  profile 

#  DBL:  position  of  beam  waist 

#  DBL:  divergence  of  beam 

#  DBL:  duration  of  one  pulse  [s] 

#  DBL:  period  of  pulse  train  [s] 

#  DBL:  time  that  emitter  actually  comes  on  [s] 

#  DBL:  time  that  emitter  turns  off  [s] 

#  DBL:  minimum  timestep  when  emitter  is  on  [s] 

#  DBL:  minimum  timestep  when  emitter  is  off  [s] 

#  DBL:  timestep  stretch  ratio  when  emitter  is  on 

#  DBL:  timestep  stretch  ratio  when  emitter  is  off 

#  STR:  name  of  sar  file  to  use 

#  STR:  name  of  file  with  temporal  pulse  profile 

#  STR:  name  of  power  spectrum  (power  vs  wavelength)  file 

#  STR:  name  of  beam  profile  file 

#  STR:  name  of  focus  vs  time  file 


Figure  12.2:  Example  of  the  emitter  configuration  file 


#KeyValue  Standard  Emitter  configuration 


#  BaseEmitter 

EmitterType 

PulseType 

SpectrumType 

Prof ileType 

PeakPower 

MinWaveLength 

rce  [1/cm] 

MaxWaveLength 

rce  [1/cm] 


Parameters 


=  "STING1 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 


=  VALUE 


#  STR:  emitter  type  (Scattering,  Propagation) 

#  INT:  pulse  type  (1  Single,  2  Multiple,  3  USER) 

#  INT:  for  broad  band  source,  type  of  wavelength  spectrum  (1  Black  Body,  2  USER) 

#  INT:  profile  type  (1  Gaussian,  2  Tophat,  3  Annular,  4  USER) 

#  DBL:  peak  power  of  emitter  [W] 

#  DBL:  minum  wavelength  for  broad  band  source,  taken  to  be  the  wavelength  of  laser  sou 

#  DBL:  minum  wavelength  for  broad  band  source,  taken  to  be  the  wavelength  of  laser  sou 


FocusType  =  VALUE 
BeamDiameter  =  VALUE 
BeamWaistPosition  =  VALUE 


TimeStepType 
dtMinOn 
dtMinOf f 
StretchOn 
StretchOf f 


=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 


PulseDuration  =  VALUE 

PulsePeriod  =  VALUE 

1  Hz,  10%  duty  cyle  train 
StartTime  =  VALUE 

StopTime  =  VALUE 


#  ScatteringEmitter  Parameters 
PhotonPacketDensity  =  VALUE 
Seed  =  VALUE 

NoiseReductionFactor  =  VALUE 
density  calculation 


#  BPEmitter  Parameters 

#NumPropSteps  =  VALUE 

zResFactor  =  VALUE 

rResFactor  =  VALUE 

OpticsFile  =  "STRING" 

#  HankelEmitter  Parameters 

p  =  VALUE 

numBesselFunctions  =  VALUE 

BesselPrecision  =  "STRING" 


#  INT:  focus  type  (1  Fixed,  2  Moving  [thermal  lens]) 

#  DBL:  beam  diameter  [cm]  Note:  this  is  1/e  for  gaussian  beam  profile 

#  DBL:  position  of  beam  waist 

#  INT:  time  step  type  (1  Fixed,  2  Adaptive) 

#  DBL:  minimum  timestep  when  emitter  is  on  [s] 

#  DBL:  minimum  timestep  when  emitter  is  off  [s] 

#  DBL:  timestep  stretch  ratio  when  emitter  is  on 

#  DBL:  timestep  stretch  ratio  when  emitter  is  off 

#  DBL:  duration  of  one  pulse  [s] 

#  DBL:  period  of  pulse  train  [s]  so,  if  period  is  1,  and  duration  is  .1,  this  would  be 

#  DBL:  time  that  emitter  actually  comes  on  [s] 

#  DBL:  time  that  emitter  turns  off  [s] 


#  DBL:  number  of  photons  per  square  cm  to  launch  in  beam  profile 

#  INT:  (OPTIONAL)  default:  1;  seed  value  to  initialize  random  number  generator 

#  INT:  (OPTIONAL)  default:  0;  the  number  of  interier  radial  rings  to  average  in  power 


#  INT:  number  z  steps  for  propagator  to  take  between  thermal  z-slices  (DEPRECIATED) 

#  INT :  replaces  NumPropSteps 

#  INT:  same  as  zResFactor,  but  for  r  direction 

#  STR:  Gaussian  Beam  Propagation  Optics  Configuration  File  (could  just  specify  Eye) 

#  INT:  order  of  Bessel  function  to  use  as  propagation  basis 

#  INT :  number  of  pth-order  Bessel  functions  to  use  in  basis 

#  STR:  precision  to  use  in  bessel  function  evaluations  (single  or  double 


Figure  12.3:  Example  of  the  standard  emitter  configuration  file 
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#KeyValue  layer  configuration  file 


LayerType 

Description 

Thickness 

Density 

Specif icHeat 

Conductivity 

ConvHeatTransRate 

Emissivity 

BloodFlowRate 


VALUE 

"STRING' 

VALUE 

VALUE 

VALUE 

VALUE 

VALUE 

VALUE 

VALUE 


INT: 
STR: 
DBL : 
DBL : 
DBL: 
DBL: 
DBL: 
DBL: 
DBL: 


layer  type  (see  layertypes.txt  for  list) 

layer  description  (retinalEP,  epidermis,  etc) 

layer  thickness  [cm] 

layer  density  [g/cmA3] 

specific  heat  of  layer  [J/g/degC] 

conductivity  [ J/s/cm/degC] 

convective  heat  transfer  rate  [ J/s/cm2/degc] 
emissivity  of  tissue  (0  -  1) 
blood  flow  rate  [g/ (cmA3*s) ] 


#Refractive  index  and  dn/dt  vs  wavelength  -  STR  ("DBL  DBL  DBL"):  — >  "\lambda  n  dn/dt" 
Refractivelndex [0]  =  "VALUE  VALUE  VALUE" 

Refractivelndex [1]  =  "VALUE  VALUE  VALUE" 


#Absorption  Coefficient  vs  wavelength 
Absorption [0]  =  "VALUE  VALUE" 

Absorption [ 1 ]  =  "VALUE  VALUE" 

#Scattering  Anisotropy  vs  wavelength 
Anisotropy [0]  =  "VALUE  VALUE" 

Anisotropy [ 1 ]  =  "VALUE  VALUE" 

♦Scattering  Coefficients  vs  wavelengt 
Scattering [0]  =  "VALUE  VALUE" 

Scattering [1]  =  "VALUE  VALUE" 


-  STR  ("DBL  DBL"):  — >  "\lambda  \mu_{a) "  [1/cm] 


-  STR  ("DBL  DBL"):  scattering  angle  cos  by  wavelength  "\lambda  \cos(\theta) 


-  STR  ("DBL  DBL"):  scattering  coef f iecients  by  wavelength  " \lambda  \mu { s ) " 


♦Reflectance  values  vs  wavelentgh 
Reflectance [0]  =  "VALUE 

Ref lectance [ 1 ]  =  "VALUE 


VALUE " 
VALUE " 


-  STR  ("DBL  DBL"):  0-1 


Temp [ 0 ] 
A  [0  ] 
Ea  [0] 
Temp [ 1 ] 
A  [1  ] 
Ea  [1] 


=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 
=  VALUE 


♦  DBL:  temperature  for  rate  ceofficients 

♦  DBL:  A  rate  coef.  [1/s] 

♦  DBL:  Ea  rate  coef.  [J/mole] 

♦  DBL:  temperature  for  rate  ceofficients 

♦  DBL:  A  rate  coef.  [1/s] 

♦  DBL:  Ea  rate  coef.  [J/mole] 


Figure  12.4:  Example  of  the  layer  configuration  file 


86 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


Acknowledgements 


This  work  was  sponsored  by  the  USAF  Research  Laboratory  and  the  Air  Force  Office  of  Scientific  Research.  C.D. 
Clark,  L.  J.  Irvin  and  G.D.  Buffington  acknowledge  the  support  of  USAF  Research  Laboratory  Contract  F44624- 
02-D7003.  RD.S.Maseberg  and  J.  Stolarski  acknowledge  the  support  of  the  Air  Force  Research  Laboratory  Human 
Effectiveness  Directorate  Consortium  Research  Fellows  Program.  The  ideas  and  opinions  presented  here  are  those  of 
the  authors  and  not  those  of  the  US  Air  Force  or  the  Department  of  Defense.  The  authors  wish  to  thank  Bo  Chen  of 
the  University  of  Texas  at  Austin  for  useful  discussions  in  the  refinement  of  thermal  models  and  boundary  condition 
parameters,  and  for  providing  validation  data. 


87 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


88 

DISTRIBUTION  A,  APPROVED  FOR  PUBLIC  RELEASE 


References 


[1]  Ronald  G.  Hadley.  Wide-angle  beam  propagation  using  pade  approximant  operators.  OPTICS  LETTERS, 
17(20):  1426-1428,  1992. 

[2]  Sami  T.  Hendow  and  Sami  A.  Shakir.  Recursive  numerical  solution  for  nonlinear  wave  propagation  in  fibers  and 
cylindrically  symmetric  systems.  Applied  Optics,  25(1 1):1759-1764,  1986. 

[3]  M.  A.  Mainster,  T.  J.  White,  J.  H.  Tips,  and  P.  W.  Wilson.  Transient  thermal  behavior  in  biological  systems. 
Bulletin  of  Math  Biophysics,  32:303-314,  1970. 

[4]  Peter  W.  Milonni  and  Joseph  H.  Eberly.  Lasers.  Wiley-Interscience,  first  edition,  1998. 

[5]  M.  Necati  Ozisik.  Boundary  value  problems  of  heat  conduction.  International  Textbook  Company,  1968. 

[6]  M.  Necati  Ozisik.  Heat  Conduction.  John  Wiley  &  Sons,  1980. 

[7]  M.  Necati  Ozisik.  Finite  Difference  Methods  in  Heat  Transfer.  CRC  Press,  1994. 

[8]  D.  W.  Peaceman  and  H.  H.  Rachford.  The  numerical  solution  of  parabolic  and  elliptic  differential  equations. 
Journal  of  the  Society  for  Industrial  and  Applied  Mathematics,  3:28-41,  1955. 

[9]  W.  H.  Press,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P.  Flannery.  Numerical  Recipes  in  C++,  The  Art  of 
Scientific  Computing.  Cambridge  University  Press,  New  York,  second  edition,  2002. 

[10]  G.  D.  Smith.  Numerical  Solution  of  Partial  Differential  Equations:  Finite  Difference  Methods.  Oxford  Applied 
Mathematics  and  Computing  Science  Series.  Oxford  University  Press,  Oxford,  thrid  edition,  1998. 

[11]  A.  Taflove  and  S.  C.  Hagness.  Computational  Electrodynamics  -  The  Finite-Difference  Time-Domain  method. 
Artech  House,  Norwood,  MA,  second  edition,  2000. 

[12]  Lihong  Wang,  Steven  L.  Jacques,  and  Liqiong  Zheng.  Mcml  -  monte  carlo  modeling  of  light  transport  in  multi¬ 
layered  tissues.  Computer  Methods  and  Programs  in  Biomedicine,  47:131-146,  1995. 

[13]  A.  J.  Welch  and  M.  J.  C.  van  Gemert.  Optical-Thermal  Response  of  Laser-Irradiated  Tissue.  Lasers,  Photonics, 
and  Electro-Optics.  Plenum  Press,  New  York,  first  edition,  1995. 

[14]  K.  S.  Yee.  Numerical  solution  of  initial  boundary  value  problems  involving  maxwell’s  equations  in  isotropic 
media.  IEEE  Transactions  Antennas  and  Propagation,  14:302-307,  1966. 


89 

DISTRIBUTION  A.  APPROVED  FOR  PUBLIC  RELEASE 


